New config parameter: geom_spread_model

CHANGELOG.md

@@ -18,6 +18,7 @@ Earthquake source parameters from P- or S-wave displacement spectra
     `noise_pre_time` and `signal_pre_time`, respectively (see pull request [#9])
   - Config parameter `rps_from_focal_mechanism` renamed to
     `rp_from_focal_mechanism` (see pull request [#9])
+  - New config parameter: `geom_spread_model` (see issue [#8])
   - Config parameters `PLOT_SHOW`, `PLOT_SAVE` and `PLOT_SAVE_FORMAT` are now
     lowercase (`plot_show`, `plot_save` and `plot_save_format`)
   - New, optional, general config parameters for specyfing author and agency
@@ -34,6 +35,11 @@ Earthquake source parameters from P- or S-wave displacement spectra
   - It is now possible to provide different vp and vs velocities, close to the
     source and close to the stations (see the new config options above and
     issue [#5])
+  - Possibility to choose a geometrical spreading model between
+    (see issue [#8]):
+    - rⁿ (default: n=1 – body waves)
+    - Boatwright et al. (2002): "r" below a cutoff distance, frequency-dependent
+      above the cutoff distance
   - QuakeML output (when using QuakeML input)
   - Removed option to read event information and traces from a pickle file
     (rarely used)
@@ -325,5 +331,6 @@ Initial Python port.
 [#3]: https://github.com/SeismicSource/sourcespec/issues/3
 [#5]: https://github.com/SeismicSource/sourcespec/issues/5
 [#6]: https://github.com/SeismicSource/sourcespec/issues/6
+[#8]: https://github.com/SeismicSource/sourcespec/issues/8
 [#9]: https://github.com/SeismicSource/sourcespec/issues/9
 [#10]: https://github.com/SeismicSource/sourcespec/issues/10

doc/source_spec.rst

@@ -28,7 +28,7 @@ propagation term (geometric and anelastic attenuation of body waves):
 .. math::
 
    S(f) =
-          \frac{1}{r}
+          \frac{1}{\mathcal{G}(r)}
           \times
           \frac{2 R_{\Theta\Phi}}
                {4 \pi \rho_h^{1/2} \rho_r^{1/2} c_h^{5/2} c_r^{1/2}}
@@ -41,8 +41,8 @@ propagation term (geometric and anelastic attenuation of body waves):
 
 where:
 
-- :math:`r` is the hypocentral distance (:math:`\frac{1}{r}` is geometric
-  attenuation);
+- :math:`\mathcal{G}(r)` is the geometrical spreading coefficient (see below)
+  and :math:`r` is the hypocentral distance;
 - :math:`R_{\Theta\Phi}` is the radiation pattern coefficient for P- or S-waves
   (average or computed from focal mechanism, if available);
 - :math:`\rho_h` and :math:`\rho_r` are the medium densities at the hypocenter
@@ -55,17 +55,62 @@ where:
 - :math:`t^*` is an attenuation parameter which includes anelastic path
   attenuation (quality factor) and station-specific effects.
 
+Geometrical spreading
+---------------------
+The geometrical spreading coefficient :math:`\mathcal{G}(r)` can be defined in
+one of the following ways:
 
-In ``source_spec``, the observed spectra :math:`S(f)` are converted into
-moment magnitude :math:`M_w`.
+- :math:`\mathcal{G}(r) = r^n`: :math:`n` can be any positive number.
+  :math:`n=1` (default value) is the theoretical value for a body wave in a
+  homogeneous full-space;
+  :math:`n=0.5` is the theoretical value for a surface wave in a homogeneous
+  half-space.
 
-The first step is to multiply the spectrum for the hypocentral distance
-and convert them to seismic moment units:
+- Follwing Boatwright et al. (2002), eq. 8:
+
+  - body wave spreading (:math:`\mathcal{G}(r) = r`) for hypocentral distances
+    below a cutoff distance :math:`r_0`;
+  - frequency dependent spreading for hypocentral distances above the
+    cutoff distance :math:`r_0`.
+
+More precisely, the expression derived from Boatwright et al. (2002) is:
+
+.. math::
+
+  \mathcal{G}(r) =
+  \begin{cases}
+    r  &  r \le r_0\\
+    r_0 (r/r_0)^{\gamma (f)}  &  r > r_0
+  \end{cases}
+
+with
+
+.. math::
+
+  \gamma (f) =
+  \begin{cases}
+    0.5  &  f \le 0.20 Hz\\
+    0.5 + 2 \log (5f)  &  0.20 < f < 0.25 Hz\\
+    0.7  &  f \ge 0.25 Hz\\
+  \end{cases}
+
+Note that here we use the square root of eq. 8 in Boatwright et al. (2002),
+since we correct the spectral amplitude and not the energy.
+
+
+Building spectra
+================
+
+In ``source_spec``, the observed spectrum :math:`S(f)` is converted into
+moment magnitude units :math:`M_w`.
+
+The first step is to multiply the spectrum for the geometrical spreading
+coefficient and convert it to seismic moment units:
 
 .. math::
 
    M(f) \equiv
-   r \times
+   \mathcal{G}(r) \times
    \frac{4 \pi \rho_h^{1/2} \rho_r^{1/2} c_h^{5/2} c_r^{1/2}}
         {2 R_{\Theta\Phi}}
    \times S(f) =
@@ -75,7 +120,7 @@ and convert them to seismic moment units:
           e^{- \pi f t^*}
 
 
-Then the spectrum is converted in unities of magnitude
+Then the spectrum is converted in units of magnitude
 (the :math:`Y_{data} (f)` vector used in the inversion):
 
 .. math::
@@ -118,6 +163,10 @@ Finally coming to the following model used for the inversion:
 
 where :math:`M_w \equiv \frac{2}{3} (\log_{10} M_0 - 9.1)`.
 
+
+Inversion procedure
+===================
+
 The parameters to determine are :math:`M_w`, :math:`f_c` and :math:`t^*`.
 
 The retrieved attenuation parameter :math:`t^*` is then converted to the P- or
@@ -128,7 +177,7 @@ S-wave quality factor :math:`Q_0^{[P|S]}` using the following expression:
    Q_0^{[P|S]} = \frac{tt_{[P|S]}(r)}{t^*}
 
 where :math:`tt_{[P|S]}(r)` is the P- or S-wave travel time from source to
-station.
+station and :math:`r` is the hypocentral distance.
 
 Station-specific effects can be determined by running ``source_spec`` on several
 events and computing the average of station residuals between observed and

sourcespec/configspec.conf

@@ -226,6 +226,28 @@ rpp = float(min=0, default=0.52)
 rps = float(min=0, default=0.62)
 # Radiation pattern from focal mechanism, if available
 rp_from_focal_mechanism = boolean(default=False)
+# Geometrical spreading correction of wave amplitude.
+# Spectra will be mutliplied by this value to correct for the lost amplitude.
+# Possible options are:
+#    '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).
+#    'boatwright': "r" (body waves) geometrical spreading for hypocentral
+#                  distances below a cutoff distance; frequency-dependent
+#                  geometrical spreading for hypocentral distances above the
+#                  cutoff distance. From eq. 8 in Boatwright et al. (2002)
+#                  (except that we take here the square root of eq. 8, since we
+#                  correct the amplitude and not the energy).
+#                  You must provide the cutoff distance (see the parameter
+#                  "geom_spread_cutoff_distance" below).
+geom_spread_model = option('r_power_n', 'boatwright', default='r_power_n')
+# Exponent "n" for the "r_power_n" geometrical spreading coefficient.
+# It can be any positive float. "n" is equal to "1" (default value) for a body
+# wave in a homogeneous full-space; it is "0.5" for a surface wave in a
+# homogeneous half-space.
+geom_spread_n_exponent = float(min=0, default=1)
+# Geometrical spreading cutoff distance, in km, for the "boatwright" model:
+geom_spread_cutoff_distance = float(min=0, default=100)
 # -------- SPECTRAL MODEL PARAMETERS
 
 

sourcespec/ssp_build_spectra.py

@@ -198,6 +198,47 @@ def _check_noise_level(trace_signal, trace_noise):
         raise RuntimeError(msg)
 
 
+def _geom_spread_r_power_n(hypo_dist_in_km, exponent):
+    """rⁿ geometrical spreading coefficient."""
+    dist = hypo_dist_in_km * 1e3
+    coeff = dist**exponent
+    return coeff
+
+
+def _geom_spread_boatwright(hypo_dist_in_km, cutoff_dist_in_km, freqs):
+    """"
+    Geometrical spreading coefficient from Boatwright et al. (2002), eq. 8.
+
+    Except that we take the square root of eq. 8, since we correct amplitude
+    and not energy.
+    """
+    dist = hypo_dist_in_km * 1e3
+    cutoff_dist = cutoff_dist_in_km * 1e3
+    if dist <= cutoff_dist:
+        coeff = dist
+    else:
+        exponent = np.ones_like(freqs)
+        low_freq = freqs <= 0.2
+        mid_freq = np.logical_and(freqs > 0.2, freqs <= 0.25)
+        high_freq = freqs >= 0.25
+        exponent[low_freq] = 0.5
+        exponent[mid_freq] = 0.5 + 2*np.log(5*freqs[mid_freq])
+        exponent[high_freq] = 0.7
+        coeff = cutoff_dist * (dist/cutoff_dist)**exponent
+    return coeff
+
+
+def _geometrical_spreading_coefficient(config, spec):
+    hypo_dist_in_km = spec.stats.hypo_dist
+    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)
+    elif config.geom_spread_model == 'boatwright':
+        cutoff_dist_in_km = config.geom_spread_cutoff_distance
+        return _geom_spread_boatwright(
+            hypo_dist_in_km, cutoff_dist_in_km, spec.get_freq())
+
+
 # store log messages to avoid duplicates
 velocity_log_messages = []
 
@@ -268,11 +309,8 @@ def _smooth_spectrum(spec, smooth_width_decades=0.2):
 
 
 def _build_spectrum(config, trace):
-    stats = trace.stats
-    # normalization for the hypocentral distance
-    trace.data *= trace.stats.hypo_dist * 1000
-    # calculate fft
     spec = spectrum.do_spectrum(trace)
+    stats = trace.stats
     spec.stats.instrtype = stats.instrtype
     spec.stats.coords = stats.coords
     spec.stats.hypo = stats.hypo
@@ -286,6 +324,9 @@ def _build_spectrum(config, trace):
         _frequency_integrate(config, spec)
     # cut the spectrum
     spec = _cut_spectrum(config, spec)
+    # correct geometrical spreading
+    geom_spread = _geometrical_spreading_coefficient(config, spec)
+    spec.data *= geom_spread
     # convert to seismic moment
     coeff = _displacement_to_moment(stats, config)
     spec.data *= coeff
@@ -426,9 +467,8 @@ def build_spectra(config, st):
     Build spectra and the spec_st object.
 
     Computes P- or S-wave (displacement) spectra from
-    accelerometers and velocimeters, uncorrected for attenuation,
-    corrected for instrumental constants, normalized by
-    hypocentral distance.
+    accelerometers and velocimeters, uncorrected for anelastic attenuation,
+    corrected for instrumental constants, normalized by geometrical spreading.
     """
     logger.info('Building spectra...')
     spec_st = Stream()

sourcespec/ssp_setup.py

@@ -393,6 +393,8 @@ def _check_mandatory_config_params(config_obj):
         'rho',
         'rpp',
         'rps',
+        'geom_spread_n_exponent',
+        'geom_spread_cutoff_distance',
         'f_weight',
         'weight',
         't_star_0',