P wave type (version 2) (#9)

* Added 'rpp' configuration parameter for P-wave average radiation pattern coefficient.

* Added support for P wave type in get_radiation_coefficient function.

* Added support for P wave type in _add_hypo_dist_and_arrivals, _merge_stream and _check_sn_ratio functions.

* Added support for P wave type in _plot_trace function.

* Added support for P wave type in _compute_h, _displacement_to_moment functions.
Added support for P wave type in _check_data_len and _cut_signal_noise functions.

* Added support for P wave type in _spec_inversion function.

* Added support for P wave type in _Q0_to_t_star function.

* Simplified distinction between P and S in _displacement_to_moment function.
Write hypocentral and station velocities to log in _displacement_to_moment function.

* Add `P` option to `wave_type` in config file

* Store `vp` in hypo object

* `pre_p_time`-->`noise_pre_time`;`pre_s_time`-->`signal_pre_time`

* Remove references to "S-wave" only

   Replace them with "P- or S-wave" or with "signal".

* Improve velocity log messages

   Add statid; avoid duplicated messages.

* `rps_from_focal_mechanism`-->`rp_from_focal_mechanism`

* CHANGELOG for #9

* Update doc for P-wave inversion

Co-authored-by: Claudio Satriano <[email protected]>

Author: krisvanneste

CHANGELOG.md

@@ -1,6 +1,6 @@
 # SourceSpec
 
-Earthquake source parameters from S-wave displacement spectra
+Earthquake source parameters from P- or S-wave displacement spectra
 
 (c) 2011-2022 Claudio Satriano <[email protected]>
 
@@ -14,6 +14,10 @@ Earthquake source parameters from S-wave displacement spectra
   - Config options `vp` and `vs` have been renamed to `vp_source`
     and `vs_source` (see issue [#5])
   - New, optional, config options `vp_stations` and `vs_stations`
+  - Config options `pre_p_time` and `pre_s_time` have been renamed to
+    `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])
   - 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
@@ -24,6 +28,7 @@ Earthquake source parameters from S-wave displacement spectra
     section called `TRACE AND METADATA PARAMETERS`
 
 - Features:
+  - Support for P-wave spectral inversion (see pull request [#9])
   - 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])
@@ -310,3 +315,4 @@ 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
+[#9]: https://github.com/SeismicSource/sourcespec/issues/9

CITATION.cff

@@ -4,6 +4,6 @@ authors:
 - family-names: "Satriano"
   given-names: "Claudio"
   orcid: "https://orcid.org/0000-0002-3039-2530"
-title: "SourceSpec – Earthquake source parameters from S-wave displacement spectra"
+title: "SourceSpec – Earthquake source parameters from P- or S-wave displacement spectra"
 doi: 10.5281/ZENODO.3688587
 url: "https://sourcespec.seismicsource.org"

README.md

@@ -2,7 +2,7 @@
 
 # SourceSpec
 
-Earthquake source parameters from S-wave displacement spectra
+Earthquake source parameters from P- or S-wave displacement spectra
 
 [![PyPI-badge]][PyPI-link]
 [![license-badge]][license-link]
@@ -97,8 +97,9 @@ repository.
 
 If you used SourceSpec for a scientific paper, please cite it as:
 
-> Satriano, C. (2022). SourceSpec – Earthquake source parameters from S-wave
-> displacement spectra (X.Y). https://doi.org/10.5281/ZENODO.3688587
+> Satriano, C. (2022). SourceSpec – Earthquake source parameters from
+> P- or S-wave displacement spectra (X.Y).
+> https://doi.org/10.5281/ZENODO.3688587
 
 Please replace `X.Y` with the SourceSpec version number you used.
 

doc/index.rst

@@ -9,17 +9,17 @@ SourceSpec documentation
 SourceSpec is a collection of command line programs (written in Python) to
 determine earthquake source parameters (seismic moment :math:`M_0`, corner
 frequency :math:`f_c`) and the inelastic attenuation term (:math:`t^*`), from
-the modeling of waveform spectra.
+the modeling of waveform displacement spectra.
 
 Other parameters (source radius :math:`r_0`, stress drop :math:`\Delta \sigma`)
-are computed from the inverted ones. The quality factor :math:`Q` is determined
-from :math:`t^*`.
+are computed from the inverted ones. The quality factor :math:`Q_0` is
+determined from :math:`t^*`.
 
 As a bonus, local magnitude :math:`M_l` is computed as well.
 
 SourceSpec is composed of the following programs:
 
-* ``source_spec``: inverts the S-wave displacement spectra from station
+* ``source_spec``: inverts the P- or S-wave displacement spectra from station
   recordings of a single event.
 * ``ssp_residuals``: computes station residuals from ``source_spec`` output.
 * ``source_model``: direct spectral modelling.

doc/source_spec.rst

@@ -4,7 +4,7 @@
 SourceSpec
 ###########
 
-Earthquake source parameters from inversion of S-wave spectra.
+Earthquake source parameters from inversion of P- or S-wave spectra.
 
 :copyright:
     2011-2022 Claudio Satriano <[email protected]>
@@ -15,7 +15,7 @@ Earthquake source parameters from inversion of S-wave spectra.
 Overview
 ========
 
-``source_spec`` inverts the S-wave displacement spectra from
+``source_spec`` inverts the P- or S-wave displacement spectra from
 station recordings of a single event.
 
 Spectral model
@@ -31,7 +31,7 @@ propagation term (geometric and anelastic attenuation of body waves):
           \frac{1}{r}
           \times
           \frac{2 R_{\Theta\Phi}}
-               {4 \pi \rho_h^{1/2} \rho_r^{1/2} \beta_h^{5/2} \beta_r^{1/2}}
+               {4 \pi \rho_h^{1/2} \rho_r^{1/2} c_h^{5/2} c_r^{1/2}}
           \times
           M_O
           \times
@@ -43,11 +43,11 @@ where:
 
 - :math:`r` is the hypocentral distance (:math:`\frac{1}{r}` is geometric
   attenuation);
-- :math:`R_{\Theta\Phi}` is the radiation pattern coefficient for S-waves
+- :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
   and at the receiver, respectively;
-- :math:`\beta_h` and :math:`\beta_r` are the S-wave velocities at the hypocenter
+- :math:`c_h` and :math:`c_r` are the P- or S-wave velocities at the hypocenter
   and at the receiver, respectively;
 - :math:`M_O` is the seismic moment;
 - :math:`f` is the frequency;
@@ -66,7 +66,7 @@ and convert them to seismic moment units:
 
    M(f) \equiv
    r \times
-   \frac{4 \pi \rho_h^{1/2} \rho_r^{1/2} \beta_h^{5/2} \beta_r^{1/2}}
+   \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) =
           M_O \times
@@ -120,14 +120,15 @@ where :math:`M_w \equiv \frac{2}{3} (\log_{10} M_0 - 9.1)`.
 
 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 S-wave
-quality factor :math:`Q_0` using the following expression:
+The retrieved attenuation parameter :math:`t^*` is then converted to the P- or
+S-wave quality factor :math:`Q_0^{[P|S]}` using the following expression:
 
 .. math::
 
-   Q_0 = \frac{tt_S(r)}{t^*}
+   Q_0^{[P|S]} = \frac{tt_{[P|S]}(r)}{t^*}
 
-where :math:`tt_S(r)` is the S-wave travel time from source to station.
+where :math:`tt_{[P|S]}(r)` is the P- or S-wave travel time from source to
+station.
 
 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

@@ -103,22 +103,23 @@ p_arrival_tolerance = float(default=4.0)
 s_arrival_tolerance = float(default=4.0)
 
 # Start time (in seconds) of the noise window, respect to the P arrival time
-pre_p_time = float(default=6.0)
+noise_pre_time = float(default=6.0)
 
-# Start time (in seconds) of the S-wave window, respect to the S arrival time
-pre_s_time = float(default=1.0)
+# Start time (in seconds) of the signal window, respect to the P or S arrival
+# times (see "wave_type" below)
+signal_pre_time = float(default=1.0)
 
-# Length (in seconds) for both noise and S-wave windows
+# Length (in seconds) for both noise and signal windows
 win_length = float(default=5.0)
 # -------- TIME WINDOW PARAMETERS
 
 
 # SPECTRUM PARAMETERS --------
-# Wave type to analyse: 'S', 'SH' or 'SV'
+# Wave type to analyse: 'P', 'S', 'SH' or 'SV'
 # If 'SH' or 'SV' are selected, traces are rotated in the radial-transverse
 # system. Transverse component is used for 'SH', radial (and vertical)
 # components are used for 'SV'
-wave_type = option('S', 'SH', 'SV', default='S')
+wave_type = option('P', 'S', 'SH', 'SV', default='S')
 
 # Integrate in time domain (default: integration in spectral domain)
 time_domain_int = boolean(default=False)
@@ -134,7 +135,7 @@ taper_halfwidth = float(0, 0.5, default=0.05)
 # this window length, so that the spectral
 # sampling is fixed to 1/spectral_win_length.
 # Comment out (or set to None) to use
-# S-wave window as spectral window length.
+# signal window as spectral window length.
 spectral_win_length = float(default=None)
 
 # Spectral smoothing window width in frequency decades
@@ -217,10 +218,12 @@ vs_stations = float(default=None)
 NLL_model_dir = string(default=None)
 # Density (kg/m3):
 rho = float(default=2500)
+# P-wave average radiation pattern coefficient:
+rpp = float(default=0.52)
 # S-wave average radiation pattern coefficient:
 rps = float(default=0.62)
 # Radiation pattern from focal mechanism, if available
-rps_from_focal_mechanism = boolean(default=False)
+rp_from_focal_mechanism = boolean(default=False)
 
 # Weighting type: 'noise', 'frequency' or 'no_weight'
 weighting = option('noise', 'frequency', 'no_weight', default='noise')

sourcespec/source_spec.py

@@ -1,7 +1,7 @@
 # -*- coding: utf-8 -*-
 # SPDX-License-Identifier: CECILL-2.1
 """
-Earthquake source parameters from inversion of S-wave spectra.
+Earthquake source parameters from inversion of P- or S-wave spectra.
 
 :copyright:
     2012 Claudio Satriano <[email protected]>

sourcespec/ssp_build_spectra.py

@@ -102,6 +102,9 @@ def _compute_h(spec_st, code, wave_type='S'):
         # only use radial and, optionally, vertical component for SV
         if wave_type == 'SV' and channel[-1] == 'T':
             continue
+        # only use vertical component for P
+        if wave_type == 'P' and channel[-1] != 'Z':
+            continue
         if spec_h is None:
             spec_h = spec.copy()
             spec_h.data = np.power(spec_h.data, 2)
@@ -118,10 +121,13 @@ def _compute_h(spec_st, code, wave_type='S'):
 
 def _check_data_len(config, trace):
     traceId = trace.get_id()
-
     trace_cut = trace.copy()
-    t1 = trace.stats.arrivals['S1'][1]
-    t2 = trace.stats.arrivals['S2'][1]
+    if config.wave_type[0] == 'S':
+        t1 = trace.stats.arrivals['S1'][1]
+        t2 = trace.stats.arrivals['S2'][1]
+    elif config.wave_type[0] == 'P':
+        t1 = trace.stats.arrivals['P1'][1]
+        t2 = trace.stats.arrivals['P2'][1]
     trace_cut.trim(starttime=t1, endtime=t2, pad=True, fill_value=0)
     npts = len(trace_cut.data)
     if npts == 0:
@@ -130,7 +136,7 @@ def _check_data_len(config, trace):
         raise RuntimeError(msg)
     nzeros = len(np.where(trace_cut.data == 0)[0])
     if nzeros > npts/4:
-        msg = '{}: S-wave window is zero for more than 25%: skipping trace'
+        msg = '{}: Signal window is zero for more than 25%: skipping trace'
         msg = msg.format(traceId)
         raise RuntimeError(msg)
 
@@ -147,8 +153,12 @@ def _cut_signal_noise(config, trace):
         _time_integrate(config, trace_noise)
 
     # trim...
-    t1 = trace.stats.arrivals['S1'][1]
-    t2 = trace.stats.arrivals['S2'][1]
+    if config.wave_type[0] == 'S':
+        t1 = trace.stats.arrivals['S1'][1]
+        t2 = trace.stats.arrivals['S2'][1]
+    elif config.wave_type[0] == 'P':
+        t1 = trace.stats.arrivals['P1'][1]
+        t2 = trace.stats.arrivals['P2'][1]
     trace_signal.trim(starttime=t1, endtime=t2, pad=True, fill_value=0)
     # Noise time window for weighting function:
     noise_t1 = trace.stats.arrivals['N1'][1]
@@ -188,21 +198,37 @@ def _check_noise_level(trace_signal, trace_noise):
         raise RuntimeError(msg)
 
 
+# store log messages to avoid duplicates
+velocity_log_messages = []
+
+
 def _displacement_to_moment(stats, config):
     """
     Return the coefficient for converting displacement to seismic moment.
 
     From Aki&Richards,1980
     """
-    vs_hypo = config.hypo.vs
-    vs_station = get_vel(
+    phase = config.wave_type[0]
+    if phase == 'S':
+        v_hypo = config.hypo.vs
+    elif phase == 'P':
+        v_hypo = config.hypo.vp
+    v_station = get_vel(
         stats.coords.longitude, stats.coords.latitude, -stats.coords.elevation,
-        'S', config)
-    vs_hypo *= 1000.
-    vs_station *= 1000.
-    vs3 = vs_hypo**(5./2) * vs_station**(1./2)
-    rps = get_radiation_pattern_coefficient(stats, config)
-    coeff = 4 * math.pi * vs3 * config.rho / (2 * rps)
+        phase, config)
+    specid = '.'.join((
+        stats.network, stats.station, stats.location, stats.channel))
+    msg = '{}: V{}_hypo: {:.2f} km/s, V{}_station: {:.2f} km/s'.format(
+            specid, phase.lower(), v_hypo, phase.lower(), v_station)
+    global velocity_log_messages
+    if msg not in velocity_log_messages:
+        logger.info(msg)
+        velocity_log_messages.append(msg)
+    v_hypo *= 1000.
+    v_station *= 1000.
+    v3 = v_hypo**(5./2) * v_station**(1./2)
+    rp = get_radiation_pattern_coefficient(stats, config)
+    coeff = 4 * math.pi * v3 * config.rho / (2 * rp)
     return coeff
 
 
@@ -399,7 +425,7 @@ def build_spectra(config, st):
     """
     Build spectra and the spec_st object.
 
-    Computes S-wave (displacement) spectra from
+    Computes P- or S-wave (displacement) spectra from
     accelerometers and velocimeters, uncorrected for attenuation,
     corrected for instrumental constants, normalized by
     hypocentral distance.

sourcespec/ssp_data_types.py

@@ -93,7 +93,8 @@ def _check_minmax(self, minmax):
             return minmax
 
     def _Qo_to_t_star(self):
-        travel_time = self.spec.stats.travel_times['S']
+        phase = self.config.wave_type[0]
+        travel_time = self.spec.stats.travel_times[phase]
         t_star_bounds = travel_time/self.config.Qo_min_max
         return sorted(t_star_bounds)
 

sourcespec/ssp_inversion.py

@@ -274,7 +274,7 @@ def _spec_inversion(config, spec, spec_weight):
 
     # additional parameters, computed from fc, Mw and t_star
     vs = config.vs_source
-    travel_time = spec.stats.travel_times['S']
+    travel_time = spec.stats.travel_times[config.wave_type[0]]
     # seismic moment
     par['Mo'] = mag_to_moment(par['Mw'])
     # source radius in meters

sourcespec/ssp_parse_arguments.py

@@ -38,11 +38,11 @@ def _get_description(progname):
     if progname == 'source_spec':
         nargs = '+'
         description = 'Estimation of seismic source parameters from '
-        description += 'inversion of S-wave spectra.'
+        description += 'inversion of P- or S-wave displacement spectra.'
         epilog = ''
     elif progname == 'source_model':
         nargs = '*'
-        description = 'Direct modeling of S-wave spectra.'
+        description = 'Direct modeling of P- or S-wave displacement spectra.'
         epilog = 'Several values of moment magnitude, seismic moment, t-star\n'
         epilog += 'and alpha can be specified using a comma-separated list,'
         epilog += 'eg.:\n'

sourcespec/ssp_plot_traces.py

@@ -198,10 +198,14 @@ def _plot_trace(config, trace, ntraces, tmax,
         ax.add_patch(rect)
     except KeyError:
         pass
-    # S-wave window
-    S1 = trace.stats.arrivals['S1'][1] - trace.stats.starttime
-    S2 = trace.stats.arrivals['S2'][1] - trace.stats.starttime
-    rect = patches.Rectangle((S1, 0), width=S2-S1, height=1,
+    # Signal window
+    if config.wave_type[0] == 'S':
+        t1 = trace.stats.arrivals['S1'][1] - trace.stats.starttime
+        t2 = trace.stats.arrivals['S2'][1] - trace.stats.starttime
+    elif config.wave_type[0] == 'P':
+        t1 = trace.stats.arrivals['P1'][1] - trace.stats.starttime
+        t2 = trace.stats.arrivals['P2'][1] - trace.stats.starttime
+    rect = patches.Rectangle((t1, 0), width=t2-t1, height=1,
                              transform=trans, color='yellow',
                              alpha=0.5, zorder=-1)
     ax.add_patch(rect)
@@ -236,7 +240,7 @@ def _add_labels(axes, plotn, ncols):
 
 def _trim_traces(config, st):
     for trace in st:
-        t1 = (trace.stats.arrivals['P'][1] - config.pre_p_time)
+        t1 = (trace.stats.arrivals['P'][1] - config.noise_pre_time)
         t2 = (trace.stats.arrivals['S'][1] + 3 * config.win_length)
         trace.trim(starttime=t1, endtime=t2, pad=True, fill_value=0)