diff --git a/CHANGELOG.md b/CHANGELOG.md index bf54a999..df8b16fc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,6 +16,9 @@ Copyright (c) 2011-2025 Claudio Satriano ### Processing +- New option `r_power_n_segmented` for the `geom_spread_model` config parameter + to use a segmented geometrical spreading model with different powers for + different distance ranges - New config parameter `refine_theoretical_arrivals` to refine the theoretical P and S arrival times using a simple autopicker based on the smoothed envelope of the trace @@ -35,6 +38,8 @@ Copyright (c) 2011-2025 Claudio Satriano ### Config file +- New option `r_power_n_segmented` for the `geom_spread_model` config parameter +- New config parameters: `geom_spread_n_exponents`, `geom_spread_n_distances` - New config parameters: `refine_theoretical_arrivals`, `autopick_freqmin`, `autopick_debug_plot` - New config parameter `clipping_min_amplitude_ratio` diff --git a/sourcespec/config_files/configspec.conf b/sourcespec/config_files/configspec.conf index 18d4effd..4706775a 100644 --- a/sourcespec/config_files/configspec.conf +++ b/sourcespec/config_files/configspec.conf @@ -356,6 +356,12 @@ rho_stations = float(min=0, default=None) # 'r_power_n': "r" to the power of "n" (rⁿ). # You must provide the value of the exponent "n" # (see "geom_spread_n_exponent" below). +# 'r_power_n_segmented': piecewise continuous rⁿ function, +# cf. Atkinson & Boore (1995), Boore (2003) +# You must provide the values of the exponent "n" +# (see "geom_spread_n_exponents" below), as well as +# the distances separating the segments +# (see "geom_spread_n_disntaces" below) # 'boatwright': "r" (body waves) geometrical spreading for hypocentral # distances below a cutoff distance; frequency-dependent # geometrical spreading above the cutoff distance (Boatwright @@ -363,12 +369,20 @@ rho_stations = float(min=0, default=None) # "geom_spread_cutoff_distance" below). This coefficient can # be a valid choice for regional distances (up to 200 km), # where S-waves, Lg waves and surface waves are mixed. -geom_spread_model = option('r_power_n', 'boatwright', default='r_power_n') +geom_spread_model = option('r_power_n', 'r_power_n_segmented', 'boatwright', default='r_power_n') # Exponent "n" for the "r_power_n" geometrical spreading coefficient (positive # float). Examples: # geom_spread_n_exponent = 1 (default, body wave in a homogeneous full-space) # geom_spread_n_exponent = 0.5 (surface wave in a homogeneous half-space) geom_spread_n_exponent = float(min=0, default=1) +# Exponents and distances (in km) separating segments with different exponents +# for the "r_power_n_segmented" geometrical spreading function. +# The number of exponents must be one more than the number of distances. +# Example: +# geom_spread_n_exponents = [1., 0., 0.5] +# geom_spread_n_distances = [70, 130] +geom_spread_n_exponents = float_list(min=0, default=None) +geom_spread_n_distances = float_list(default=None) # Geometrical spreading cutoff hypocentral distance, in km, for the # "boatwright" model: geom_spread_cutoff_distance = float(min=0.01, default=50) diff --git a/sourcespec/ssp_build_spectra.py b/sourcespec/ssp_build_spectra.py index d84ee66d..b01af23d 100644 --- a/sourcespec/ssp_build_spectra.py +++ b/sourcespec/ssp_build_spectra.py @@ -29,7 +29,8 @@ from sourcespec.ssp_setup import ssp_exit from sourcespec.ssp_util import ( smooth, cosine_taper, moment_to_mag, MediumProperties, - geom_spread_r_power_n, geom_spread_boatwright, geom_spread_teleseismic) + geom_spread_r_power_n, geom_spread_r_power_n_segmented, + geom_spread_boatwright, geom_spread_teleseismic) from sourcespec.ssp_process_traces import filter_trace from sourcespec.ssp_correction import station_correction from sourcespec.ssp_radiation_pattern import get_radiation_pattern_coefficient @@ -273,6 +274,17 @@ def _geometrical_spreading_coefficient(config, spec): if config.geom_spread_model == 'r_power_n': exponent = config.geom_spread_n_exponent return geom_spread_r_power_n(hypo_dist_in_km, exponent) + if config.geom_spread_model == 'r_power_n_segmented': + exponents = config.geom_spread_n_exponents + hinge_distances = config.geom_spread_n_distances + if len(hinge_distances) != len(exponents) - 1: + raise ValueError( + f'The number of exponents must be one more than the number of ' + f'hinge distances. You provided {len(exponents)} exponents and ' + f'{len(hinge_distances)} hinge distances' + ) + return geom_spread_r_power_n_segmented(hypo_dist_in_km, exponents, + hinge_distances) if config.geom_spread_model == 'boatwright': cutoff_dist_in_km = config.geom_spread_cutoff_distance return geom_spread_boatwright( diff --git a/sourcespec/ssp_util.py b/sourcespec/ssp_util.py index 2e6eb670..6e4f9cf7 100644 --- a/sourcespec/ssp_util.py +++ b/sourcespec/ssp_util.py @@ -267,6 +267,47 @@ def geom_spread_r_power_n(hypo_dist_in_km, exponent): return dist**exponent +def geom_spread_r_power_n_segmented(hypo_dist_in_km, exponents, + hinge_distances): + """ + Geometrical spreading function defined as piecewise continuous powerlaw, + as defined in Boore (2003), eq. 9 + + :param hypo_dist_in_km: Hypocentral distance (km). + :type hypo_dist_in_km: float + :param exponents: Exponents for different powerlaw segments + :type hinge_distances: numpy.ndarray + :param hinge_distances: Distances between different powerlaw segments + :type exponent: numpy.ndarray + :return: Geometrical spreading correction (for distance in m) + :rtype: float + """ + if np.isscalar(hypo_dist_in_km): + hypo_dist_in_km = np.array([hypo_dist_in_km], dtype='float') + is_scalar = True + else: + hypo_dist_in_km = np.asarray(hypo_dist_in_km, dtype='float') + is_scalar = False + Rref = 1. + hinge_distances = np.hstack([[Rref], hinge_distances or []]) + exponents = -np.asarray(exponents) + # Do not allow distances less than Rref (1 km) + hypo_dist_in_km = np.maximum(Rref, hypo_dist_in_km) + Zhinges = (hinge_distances[:-1] / hinge_distances[1:]) ** exponents[:-1] + Zhinges = np.cumprod(Zhinges) + R0, p0 = hinge_distances[0], exponents[0] + Z = (R0 / hypo_dist_in_km) ** p0 + for n in range(1, len(hinge_distances)): + Rn, pn = hinge_distances[n], exponents[n] + idxs = hypo_dist_in_km > Rn + Z[idxs] = Zhinges[n-1] * ((Rn / hypo_dist_in_km[idxs]) ** pn) + # Convert spreading correction to metric distance + Z *= 1e3 + if is_scalar: + Z = Z[0] + return Z + + def _boatwright_above_cutoff_dist(freqs, cutoff_dist, dist): """ Geometrical spreading coefficient from Boatwright et al. (2002), eq. 8.