Merge pull request #62 from krisvanneste/segmented_r_power_n

Segmented r_power_n

CHANGELOG.md

@@ -16,6 +16,9 @@ Copyright (c) 2011-2025 Claudio Satriano <[email protected]>
 
 ### 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 <[email protected]>
 
 ### 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`

sourcespec/config_files/configspec.conf

@@ -356,19 +356,33 @@ 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
 #                   et al., 2002). You must provide the cutoff distance (see
 #                   "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)

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(

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.