compute rfs - GeoscienceAustralia/hiperseis GitHub Wiki

Compute Receiver Functions (RFs)

The RF workflow consists of:

  • Computing RFs
  • Quality labelling and reporting

In addition to documentation below, see module level README file.

Generating RFs from event seismograms

RFs are computed by taking event seismograms and job configuration parameters and running them through the generate_rf.py script:

This produces a new HDF5 file containing 3-component RFs (ZRT or LQT, depending on the user setting). This RF HDF5 file is then ingested into subsequent workflows below.

The generate_rf.py script takes the following options. See generate_rf.py --help for correct ordering of command line arguments:

  • Name of job configuration file.
  • Name of the input event waveform file.
  • Name of the output file to write RFs to.

All the details of how to compute the RFs are contained dictionary format in the JSON configuration file. At the time of this writing, the configuration fields are as follows. Note that this is not valid JSON, and not all fields are required, this just shows the key names and indicative types. See the source code documentation for full explanations of fields. Users should study some examples before trying to compose these configuration files:

{
  "filtering":
  {
    "resample_rate": float,
    "taper_limit": float,
    "filter_band": [float, float],
    "channel_pattern": str,
    "baz_range": [float, float]
  },

  "processing":
  {
    "custom_preproc":
    {
      "import": str,
      "func": str,
      "args": {}
    },
    "trim_start_time": float,
    "trim_end_time": float,
    "rotation_type": str,
    "deconv_domain": str,
    "gauss_width": float,
    "water_level": float,
    "spiking": float,
    "normalize": bool
  },

  "system":
  {
    "parallel": bool,
    "memmap": bool,
    "temp_dir": str,
    "aggressive_dispatch": bool
  }
}

Note that iterative deconvolution produces cleaner RFs, but takes much longer than time or frequency domain deconvolution.

Labelling RFs with quality metrics

Filtering RFs prior to stacking is important to obtaining useful results. However there is no consensus on what characteristics constitute 'good' quality, and it may depend on the particulars of your dataset. Therefore maximum flexibility around RF filtering is considered a desirable attribute of this workflow. Considering that some RF quality metrics can be computationally expensive to compute, the approach used here is to defer RF quality filtering to as late as possible in the processing pipeline, and to compute said quality metrics only once per RF.

Therefore computation of RF quality metrics is performed as a distinct step. This is performed by the rf_quality_filter.py script shown as the top left to right process flow in the following schematic, up to 'RFs quality labelled':

The quality labelling script simply reads all the RFs and saves a copy of them labelled with additional quality metrics in the stats field of the individual seismic traces. rf_quality_filter.py script takes the following arguments:

  • Name of input HDF5 file containing RFs generated by generate_rf.py
  • Name of output HDF5 in which to save RFs labelled with quality metrics

Report generation and H-k stacking

The file with RF quality labelling is ingested into the bulk_rf_report.py script to produce the RF report file in PDF format, as shown in the schematic above. H-k stacking is performed as an intrinsic part of the RF reporting script, since it is cheap and fast to run, and deferring computation of the H-k stack until reporting time allows for more readily customizable weighting parameters.

The bulk_rf_report.py script has various options that can be found by running bulk_rf_report.py --help. The following arguments are required for report generation:

  • Name of input HDF5 file containing RFs with quality labels, generated by rf_quality_filter.py script.
  • Name of output PDF file in which to save the report.

It is also recommended to provide at least one filtering option (unless you want to see stacks of RFs with no quality filtering) and explicit H-k stacking weightings:

  • Basic filtering found to work well with recent broadband temporal deployments: --apply-amplitude-filter
  • H-k weighting values as a triplet of numeric values, e.g. 0.34 0.33 0.33 is often used to give equal weight to all three phases included in the stacking algorithm.
⚠️ **GitHub.com Fallback** ⚠️