pyadjoint.adjoint_source_types.multitaper_misfit

Multitaper based phase and amplitude misfit and adjoint source.

authors

adjTomo Dev Team (adjtomo@gmail.com), 2022 Youyi Ruan (youyir@princeton.edu), 2016 Matthieu Lefebvre (ml15@princeton.edu), 2016 Yanhua O. Yuan (yanhuay@princeton.edu), 2015

license

BSD 3-Clause (“BSD New” or “BSD Simplified”)

Module Contents

Classes

MultitaperMisfit

A class to house the machinery of the multitaper misfit calculation. This is

Functions

calculate_adjoint_source(observed, synthetic, config, ...)

Convenience wrapper function for MTM class to match the expected format
class pyadjoint.adjoint_source_types.multitaper_misfit.MultitaperMisfit(observed, synthetic, config, windows, observed_2=None, synthetic_2=None, windows_2=None)[source]

A class to house the machinery of the multitaper misfit calculation. This is done with a class rather than a function to reduce the amount of unnecessary parameter passing between functions.

calculate_adjoint_source()[source]

Main processing function to calculate adjoint source for MTM. The logic here is that the function will perform a series of checks/ calculations to check if MTM is valid for each given window. If any check/calculation fails, it will fall back to CCTM misfit for the given window.

calculate_dd_adjoint_source()[source]

Process double difference adjoint source. Requires second set of waveforms and windows.

Note

amplitude measurement stuff has been mostly left in the function (i.e., commented out) even it is not used, so that hopefully it is easier for someone in the future to implement it if they want.

Return type

(float, np.array, np.array, dict)

Returns

(misfit_sum_p, fp, fp_2, win_stats) == ( total phase misfit for the measurement, adjoint source for first data-synthetic pair, adjoint source for second data-synthetic pair, measurement information dictionary )

calculate_mt_adjsrc(s, tapers, nfreq_min, nfreq_max, dtau_mtm, dlna_mtm, wp_w, wq_w)[source]

Calculate the adjoint source for a multitaper measurement, which tapers synthetics in various windowed frequency-dependent tapers and scales them by phase dependent travel time measurements (which incorporate the observed data).

Parameters

  • s (np.array) – synthetic data array

  • tapers (np.array) – array of DPPS windows shaped (num_taper, nlen_w)

  • nfreq_min (int) – minimum frequency for suitable MTM measurement

  • nfreq_max (int) – maximum frequency for suitable MTM measurement

  • dtau_mtm (np.array) – phase dependent travel time measurements from mtm

  • dlna_mtm (np.array) – phase dependent amplitude anomaly

  • wp_w (np.array) – phase-misfit error weighted frequency domain taper

  • wq_w (np.array) – amplitude-misfit error weighted frequency domain taper

calculate_dd_mt_adjsrc(s, s_2, tapers, nfreq_min, nfreq_max, df, dtau_mtm, dlna_mtm, wp_w, wq_w)[source]

Calculate the double difference adjoint source for multitaper measurement. Almost the same as calculate_mt_adjsrc but only addresses phase misfit and requres a second set of synthetics s_2 which is processed in the same way as the first set s

Parameters

  • s (np.array) – synthetic data array

  • s_2 (np.array) – optional 2nd set of synthetics for double difference measurements only. This will change the output

  • tapers (np.array) – array of DPPS windows shaped (num_taper, nlen_w)

  • nfreq_min (int) – minimum frequency for suitable MTM measurement

  • nfreq_max (int) – maximum frequency for suitable MTM measurement

  • df (floats) – step length of frequency bins for FFT

  • dtau_mtm (np.array) – phase dependent travel time measurements from mtm

  • dlna_mtm (np.array) – phase dependent amplitude anomaly

  • wp_w (np.array) – phase-misfit error weighted frequency domain taper

  • wq_w (np.array) – amplitude-misfit error weighted frequency domain taper

calculate_freq_domain_taper(nfreq_min, nfreq_max, df, dtau_mtm, dlna_mtm, err_dt_cc, err_dlna_cc, err_dtau_mt, err_dlna_mt)[source]

Calculates frequency domain taper weighted by misfit (either CC or MTM)

Note

Frequency-domain tapers are based on adjusted frequency band and error estimation. They are not one of the filtering processes that needs to be applied to the adjoint source but rather a frequency domain weighting function for adjoint source and misfit function.

Parameters

  • nfreq_min (int) – minimum frequency for suitable MTM measurement

  • nfreq_max (int) – maximum frequency for suitable MTM measurement

  • df (floats) – step length of frequency bins for FFT

  • dtau_mtm (np.array) – phase dependent travel time measurements from mtm

  • dlna_mtm (np.array) – phase dependent amplitude anomaly

  • err_dt_cc (float) – cross correlation time shift error

  • err_dlna_cc (float) – cross correlation amplitude anomaly error

  • err_dtau_mt (np.array) – phase-dependent timeshift error

  • err_dlna_mt (np.array) – phase-dependent amplitude error

calculate_multitaper(d, s, tapers, wvec, nfreq_min, nfreq_max, cc_tshift, cc_dlna)[source]

Measure phase-dependent time shifts and amplitude anomalies using the multitaper method

Note

Formerly mt_measure. Renamed for additional clarity and to match the CCTM function names

Parameters

  • d (np.array) – observed data array

  • s (np.array) – synthetic data array

  • tapers (np.array) – array of DPPS windows shaped (num_taper, nlen_w)

  • wvec (np.array) – angular frequency array generated from Discrete Fourier Transform sample frequencies

  • nfreq_min (int) – minimum frequency for suitable MTM measurement

  • nfreq_max (int) – maximum frequency for suitable MTM measurement

  • cc_tshift (float) – cross correlation time shift

  • cc_dlna (float) – amplitude anomaly from cross correlation

Return type

tuple of np.array

Returns

(phi_w, abs_w, dtau_w, dlna_w); (frequency dependent phase anomaly, phase dependent amplitude anomaly, phase dependent cross-correlation time shift, phase dependent cross-correlation amplitude anomaly)

calculate_mt_error(d, s, tapers, wvec, nfreq_min, nfreq_max, cc_tshift, cc_dlna, phi_mtm, abs_mtm, dtau_mtm, dlna_mtm)[source]

Calculate multitaper error with Jackknife MT estimates.

The jackknife estimator of a parameter is found by systematically leaving out each observation from a dataset and calculating the parameter estimate over the remaining observations and then aggregating these calculations.

Parameters

  • d (np.array) – observed data array

  • s (np.array) – synthetic data array

  • tapers (np.array) – array of DPPS windows shaped (num_taper, nlen_w)

  • wvec (np.array) – angular frequency array generated from Discrete Fourier Transform sample frequencies

  • nfreq_min (int) – minimum frequency for suitable MTM measurement

  • nfreq_max (int) – maximum frequency for suitable MTM measurement

  • cc_tshift (float) – cross correlation time shift

  • cc_dlna (float) – amplitude anomaly from cross correlation

  • phi_mtm (np.array) – frequency dependent phase anomaly

  • abs_mtm (np.array) – phase dependent amplitude anomaly

  • dtau_mtm (np.array) – phase dependent cross-correlation time shift

  • dlna_mtm (np.array) – phase dependent cross-correlation amplitude anomaly)

Return type

tuple of np.array

Returns

(err_phi, err_abs, err_dtau, err_dlna), (error in frequency dependent phase anomaly, error in phase dependent amplitude anomaly, error in phase dependent cross-correlation time shift, error in phase dependent cross-correlation amplitude anomaly)

calculate_freq_limits(df)[source]

Determine if a given window is suitable for multitaper measurements. If so, finds the maximum frequency range for the measurement using a spectrum of tapered synthetic waveforms

First check if the window is suitable for mtm measurements, then find the maximum frequency point for measurement using the spectrum of tapered synthetics.

Note

formerly frequency_limit. renamed to be more descriptive. also split off earlier cycle check from this function into check_sufficient_number_of_wavelengths

Parameters

df (float) – step length of frequency bins for FFT

Return type

tuple

:return (float, float, bool);

(minimumum frequency, maximum frequency, continue with MTM?)

prepare_data_for_mtm(d, tshift, dlna, window)[source]

Re-window observed data to center on the optimal time shift, and scale by amplitude anomaly to get best matching waveforms for MTM

Returns

check_time_series_acceptability(cc_tshift, nlen_w)[source]

Checking acceptability of the time series characteristics for MTM

Parameters

  • cc_tshift (float) – time shift in unit [s]

  • nlen_w (int) – window length in samples

Return type

bool

Returns

True if time series OK for MTM, False if fall back to CC

check_mtm_time_shift_acceptability(nfreq_min, nfreq_max, df, cc_tshift, dtau_mtm, sigma_dtau_mt)[source]

Check MTM time shift measurements to see if they are within allowable bounds set by the config. If any of the phases used in MTM do not meet criteria, we will fall back to CC measurement.

Note

formerly mt_measure_select, renamed for clarity

Parameters

  • nfreq_max (int) – maximum in frequency domain

  • nfreq_min (int) – minimum in frequency domain

  • df (floats) – step length of frequency bins for FFT

  • cc_tshift (float) – c.c. time shift

  • dtau_mtm (np.array) – phase dependent travel time measurements from mtm

  • sigma_dtau_mt (np.array) – phase-dependent error of multitaper measurement

Return type

bool

Returns

flag for whether any of the MTM phases failed check

rewindow(data, left_sample, right_sample, shift)[source]

Align data in a window according to a given time shift. Will not fully shift if shifted window hits bounds of the data array

Parameters

  • data (np.array) – full data array to cut with shifted window

  • left_sample (int) – left window border

  • right_sample (int) – right window border

  • shift (int) – overall time shift in units of samples

static _search_frequency_limit(is_search, index, nfreq_limit, spectra, water_threshold, c=10)[source]

Search valid frequency range of spectra. If the spectra larger than 10 * water_threshold it will trigger the search again, works like the heating thermostat.

Parameters

  • is_search (bool) – Logic switch

  • index (int) – index of spectra

  • nfreq_limit (int) – index of freqency limit searched

  • spectra (int) – spectra of signal

  • water_threshold (float) – optional triggering value to stop the search, if not given, defaults to Config value

  • c (int) – constant scaling factor for water threshold.

pyadjoint.adjoint_source_types.multitaper_misfit.calculate_adjoint_source(observed, synthetic, config, windows, observed_2=None, synthetic_2=None, windows_2=None)[source]

Convenience wrapper function for MTM class to match the expected format of Pyadjoint. Contains the logic for what to return to User.

Parameters

  • observed (obspy.core.trace.Trace) – observed waveform to calculate adjoint source

  • synthetic (obspy.core.trace.Trace) – synthetic waveform to calculate adjoint source

  • config (pyadjoint.config.ConfigCCTraveltime) – Config class with parameters to control processing

  • windows (list of tuples) – [(left, right),…] representing left and right window borders to be used to calculate misfit and adjoint sources

  • observed_2 (obspy.core.trace.Trace) – second observed waveform to calculate adjoint sources from station pairs

  • synthetic_2 (obspy.core.trace.Trace) – second synthetic waveform to calculate adjoint sources from station pairs

  • windows_2 (list of tuples) – [(left, right),…] representing left and right window borders to be tapered in units of seconds since first sample in data array. Used to window observed_2 and synthetic_2