lsforce.lsforce module
- class lsforce.lsforce.LSForce(data, data_sampling_rate, main_folder=None, method='full')[source]
Bases:
object
Class for performing force inversions.
- filter
Dictionary with keys
'freqmin'
,'freqmax'
,'zerophase'
,'periodmin'
,'periodmax'
, and'order'
specifying filter parameters- Type:
- W
Weight matrix
- Type:
2D array
- Wvec
Weight vector
- Type:
1D array
- jackknife
Dictionary with keys
'Z'
,'N'
,'E'
,'VR_all'
,'alphas'
,'num_iter'
, and'frac_delete'
containing jackknife parameters and results- Type:
- angle_magnitude
Dictionary with keys
'magnitude'
,'magnitude_upper'
,'magnitude_lower'
,'vertical_angle'
, and'horizontal_angle'
containing inversion angle and magnitude information- Type:
- G
Design matrix
- Type:
2D array
- d
Data vector
- Type:
1D array
- model
Model vector of concatenated components (n x 1) of solution
- Z
[N] Vertical force time series extracted from model (positive up)
- N
[N] North force time series extracted from model (positive north)
- E
[N] East force time series extracted from model (positive east)
- tvec
[s] Time vector for forces, referenced using zero_time (if specified)
- VR
[%] Variance reduction. Rule of thumb: This should be ~50–80%, if ~100%, solution is fitting data exactly and results are suspect. If ~5%, model may be wrong or something else may be wrong with setup
- dtorig
Original data vector
- dtnew
Modeled data vector (Gm-d)
- alpha
Regularization parameter that was used
- alphafit
Dictionary with keys
'alphas'
,'fit'
, and'size'
specifying regularization parameters tested- Type:
Create an LSForce object.
- Parameters:
data (
LSData
) – LSData object, corrected for station response but not filtereddata_sampling_rate (int or float) – [Hz] Samples per second to use in inversion. All data will be resampled to this rate, and Green’s functions will be created with this rate
main_folder (str) – If None, will use current folder
method (str) – How to parameterize the force-time function. One of ‘full’ — full inversion using Tikhonov regularization (L2 norm minimization) or ‘triangle’ — inversion parameterized using overlapping triangles, variation on method of Ekström & Stark (2013)
- forward(Z, N, E)[source]
Execute the forward problem \(\mathbf{d} = \mathbf{G}\mathbf{m}\) using a user-supplied force time series \(\mathbf{m}\) composed of components Z, N, and E.
- Parameters:
Z – [N] Vertical force time series (positive up)
N – [N] North force time series (positive north)
E – [N] East force time series (positive east)
- Returns:
[m] Stream containing synthetic data, \(\mathbf{d}\)
- Return type:
- invert(zero_time=None, impose_zero_start=False, add_to_zero=False, duration=None, jackknife=False, num_iter=200, frac_delete=0.5, alpha=None, zero_scaler=2.0, zero_start_taper_length=0, tikhonov_ratios=(1.0, 0.0, 0.0), jk_refine_alpha=False, save_matrices=False)[source]
Performs single-force inversion using Tikhonov regularization.
- Parameters:
zero_time (int or float) – [s] Optional estimated start time of real (landslide-related) part of signal, in seconds from start time of seismic data. Useful for making figures showing selected start time and also for the impose_zero_start option
impose_zero_start (bool) – Adds weighting matrix to suggest that forces tend towards zero prior to zero_time (zero_time must be defined)
add_to_zero (bool) – Adds weighting matrix to suggest that all components of force integrate to zero
duration (int or float) – Maximum duration allowed for the event, starting at zero_time if defined, otherwise starting from the beginning of the seismic data. Forces after this will tend towards zero. This helps tamp down artifacts due to edge effects, etc.
jackknife (bool) – If True, perform num_iter additional iterations of the model while randomly discarding frac_delete of the data
num_iter (int) – Number of jackknife iterations to perform
frac_delete (int or float) – Fraction (out of 1) of data to discard for each iteration, if frac_delete=1, will do leave a one out error analysis
alpha (int or float) – Set regularization parameter. If None, will search for best alpha using the L-curve method
zero_scaler (int or float) – Relative strength of zero constraint for impose_zero_start and duration options. Ranges from 0 to 10. The lower the number, the weaker the constraint. Values up to 30 are technically allowed but discouraged because high zero_scaler values risk the addition of high frequency oscillations due to the sudden release of the constraint
zero_start_taper_length (int or float) – [s] Length of taper for impose_zero_start option
tikhonov_ratios (list or tuple) – Proportion each regularization method contributes to the overall regularization effect, where values correspond to [0th order, 1st order, 2nd order]. Must sum to 1
jk_refine_alpha (bool) – Refine the alpha parameter used for each jackknife iteration by searching over order of magnitude around the best alpha for the full solution. If False, each jackknife iteration will use the same alpha as the main solution (note that this is much faster but can result in some jackknife iterations having depressed amplitudes)
save_matrices (bool) – If True, will save the inverted matrices as part of the object (Ghat, dhat, I, L1, L2) in case user wants to do additional alpha searching
- plot_fits(equal_scale=True, xlim=None)[source]
Create a plot showing the model-produced waveform fit to the data.
- plot_forces(subplots=False, xlim=None, ylim=None, same_y=True, highf_tr=None, hfshift=0.0, hfylabel=None, infra_tr=None, infra_shift=0, jackshowall=False)[source]
Plot inversion result.
- Parameters:
subplots (bool) – If True, make subplots for components, otherwise plot all on one plot
same_y (bool) – If True, use same y-axis limits for all plots
highf_tr (
Trace
) – Seismic trace with start time identical to start time of the data used in the inversionhfylabel (str) – Label used for seismic trace. If not defined, will use station name
infra_tr (
Trace
) – Infrasound trace with start time identical to start time of the data used in the inversioninfra_shift (int or float) – [s] Time shift for infrasound trace
jackshowall (bool) – If True and jackknife was run, will show all individual runs (changes subplots to True)
- Returns:
Output figure handle
- Return type:
- saverun(prefix, filepath=None, timestamp=False, figs2save=None, figs2save_names=None, save_size='light', filetype='png')[source]
Save a force inversion run for later use.
Warning
Do not expect this to work if you have the
autoreload
IPython extension enabled!- Parameters:
prefix (str) – Run name to prepend to all output files
filepath (str) – Full path to directory where all files should be saved. If None, will use self.main_folder
timestamp (bool) – Name results with current time to avoid overwriting previous results
figs2save_names (list or tuple) – Names of figures (appends to end)
save_size (str or None) – If light, does not save seismic data with object to save size. If ultralight, does not save seismic data, Green’s functions, weight matrices, design matrix, data vector, and model vector. If None, the full object is saved.
filetype (str) – Filetype given as extension, e.g. ‘png’
- setup(period_range, syngine_model=None, cps_model=None, triangle_half_width=None, source_depth=0, weights=None, noise_window_dur=None, filter_order=2, zerophase=True, skip_datafilter=False)[source]
Downloads/computes Green’s functions (GFs) and creates all matrices.
- Parameters:
syngine_model (str) – Name of Syngine model to use. If this is not None, then we calculate GFs using Syngine (preferred)
cps_model (str) – Filename of CPS model to use. If this is not None, then we calculate GFs using CPS
triangle_half_width (int or float) – [s] Half-width of triangles; only used if the triangle method is being used
weights (list or tuple or str) – If None, no weighting is applied. An array of floats with length
st_proc.count()
(and in the order of thest_proc
attribute of theLSData
object) applies manual weighting. If ‘prenoise’, uses standard deviation of a noise window defined by noise_window_dur to weight. If ‘distance’, weights by 1 / distancenoise_window_dur (int or float) – [s] Length of noise window for ‘prenoise’ weighting scheme (if not None, weights is set to ‘prenoise’)
filter_order (int) – Order of filter applied over period_range
zerophase (bool) – If True, zero-phase filtering will be used
skip_datafilter (bool) – If True, filtering will not be applied to the input data and will only be applied to the Green’s functions. This should be chosen only if the data were pre-filtered manually already with the same band as period_range so the user doesn’t want to filter them again
- lsforce.lsforce.Lcurve(fit1, size1, alphas, bestalpha=None)[source]
Plot an L-curve.
- Parameters:
fit1 (1D array) – List of residuals
size1 (1D array) – List of model norms
alphas (1D array) – List of alphas tried
bestalpha (float) – The alpha value chosen
- Returns:
figure handle
- lsforce.lsforce.find_alpha(Ghat, dhat, I, L1=0, L2=0, tikhonov_ratios=(1.0, 0.0, 0.0), rough=False, range_rough=None, int_rough=0.75, plot_Lcurve=True)[source]
Finds best regularization (trade-off) parameter alpha.
Computes model with many values of alpha, plots L-curve, and finds point of steepest curvature where slope is negative.
- Parameters:
Ghat (array) – (m x n) matrix
dhat (array) – (1 x n) array of weighted data
I (array) – Identity matrix
L1 (array) – First order roughening matrix. If 0, will use only 0th-order Tikhonov regularization
L2 (array) – Second order roughening matrix. If 0, will use only 0th-order Tikhonov regularization
tikhonov_ratios (list or tuple) – Proportion each regularization method contributes to the overall regularization effect, where values correspond to [0th order, 1st order, 2nd order]. Must sum to 1
rough (bool) – If False, will do two iterations to fine tune the alpha parameter. The second iteration searches over +/- one order of magnitude from the best alpha found from the first round. If True, time will be saved because it will only do one round of searching.
range_rough (tuple) – Lower and upper bound of range to search over in log units. If None, the program will choose a range based on the norm of
Ghat
int_rough (float) – Interval, in log units, to use for rough alpha search
plot_Lcurve (bool) – Toggle showing the L-curve plot
- Returns:
Tuple containing:
bestalpha (float) – The optimal alpha
fit1 (1D array) – List of residuals
size1 (1D array) – List of model norms
alphas (1D array) – List of alphas tried
- Return type: