Stack PRFs with Common Conversion points (CCP) method

This note will introduce how to stack PRFs with CCP method in two steps.

1. Convert PRFs from time axis to depth axis according to 1D velocity model.

2. Search and stack amplitudes in each bin with calculated pierce points falling in.

Preparation

Before the CCP stacking, we should prepare:

PRFs data with sac format

  • The SACHeader of b, delta and npts are required. please make sure that they are correct. The b denote the shift time before P arrival.

  • All of SAC files belong to a same station must be put into a same directory with name of the station name.

  • A list file with 10 columns is required in this directory, whose format as:

    • evt: the datetime part of sacfile name of each PRF.

    • phase: The phase part of sacfile name of each PRF.

    • evlat: The latitude of the event.

    • evlon: The longitude of the event.

    • evdep: The focal depth of the event.

    • dis: The epicentral distance between event and station.

    • bazi: The back-azimuth from event to station.

    • rayp: The ray parameter of the event.

    • mag: The magnitude of the event.

    • f0: The gauss factor used in computing the PRF.

A station list

A text table of station information with 3 columns:

  • station: The station name which should be the same as the directory name including PRF Files.

  • stlat: The latitude of the station.

  • stlon: The longitude of the station.

A parameter file

The parameter file include parameters used in the time-to-depth and CCP stacking. the format should follow configparser module. See Template for CCP stacking

(Optional) A lib of Ps ray parameters

If you would stack PRFs in a large depth (e.g., D410 or D660), the ray parameters of P arrival cannot represent that of Ps phase. Therefore, a library file of Ps ray parameters in different depths and epicentral distance is required in calculating Ps-P time difference. This lib file has specific binary format. You can generate it by command of gen_rayp_lib:

usage: gen_rayp_lib [-h] -d min_dis/max_dis/interval -e min_dep/max_dep/interval [-l min_layer/man_layer] [-o outpath]

Gen a ray parameter lib for Pds phases

optional arguments:
  -h, --help            show this help message and exit
  -d min_dis/max_dis/interval
                        Distance range in degree
  -e min_dep/max_dep/interval
                        Event depth range in km
  -l min_layer/man_layer
                        layers range as in km, defaults to 0/800
  -o outpath            Out path to Pds ray parameter lib

Here, we provide a ray-parameter file with distance from 30-90 degree and source depth from 0 to 600 km: psrayp.npz

(Optional) A 3D velocity model

To correct influence by velocity anomalies above the interface you wanna study, you can use a 3D velocity model to correct Ps-P time difference. The model file should be .npz format (generated by numpy.savez) with 5 fields:

  • dep: Depth axis. 1D numpy.array with shape of (len_dep,)

  • lat: Latitude axis. 1D numpy.array with shape of (len_lat,)

  • lon: Longitude axis. 1D numpy.array with shape of (len_lon,)

  • vp: P wave velocity. 3D numpy.array with shape of (len_dep, len_lat, len_lon)

  • vs: S wave velocity. 3D numpy.array with shape of (len_dep, len_lat, len_lon)

Convert time to depth

rf2depth command would convert PRFs from time axis to depth axis using above preparations:

usage: rf2depth [-h] [-d D] [-m M] [-r R] cfg_file

Convert Ps RF to depth axis

positional arguments:
  cfg_file    Path to configure file

optional arguments:
  -h, --help  show this help message and exit
  -d D        Path to 3d vel model in npz file for moveout correcting
  -m M        Folder path to 1d vel model files with staname.vel as the file name
  -r R        Path to 3d vel model in npz file for 3D ray tracing

The output structure would be saved as a .mat file, which can be read in MATLAB. Therefore, you can browse items and fields of this file using MATLAB or scipy.io.loadmat. The number of items is equal to the number of stations in Station list. Each item has 9 fields:

  • Station: The station name.

  • stalat: The Latitude of the station.

  • stalon: The Longitude of the station.

  • bazi: The back-azimuth of each event (1D array with shape of (ev_num,)).

  • rayp: The Ray-parameter of each event (1D array with shape of (ev_num,)).

  • moveout_correct: The amplitude for each PRF after time-depth conversion (2D array with shape of (ev_num, layer_num)).

  • Piercelat: Latitudes of each event at each depth (2D array with shape of (ev_num, layer_num)).

  • Piercelon: Longitudes of each event at each depth (2D array with shape of (ev_num, layer_num)).

  • StopIndex: The last layer after time-depth conversion (2D array with shape of (ev_num, layer_num)).

NOTE: The layer_num was determined by field [depth] in parameter file.

CCP Stack PRFs along a profile

The ccp_profile command provide functions to stack PRFs along a profile:

usage: ccp_profile [-h] [-t] cfg_file

Stack PRFS along a profile

positional arguments:
  cfg_file    Path to CCP configure file

optional arguments:
  -h, --help  show this help message and exit
  -t          Output as a text file

Required parameters

  • cfg_file: the parameter file introduced above. Meanwhile, some options are required:

    • FileIO section: input and output path

      • depthdat: The .mat file path output from rf2depth

      • stackfile: Path to output file after CCP stacking

      • stalist: Path to station list used to stack. If it is not exists, the stations would be limited using the width of the profile. the stations that meets the condition would be project to this profile and write to this stalist.

    • stack section: Define samples in depth

      • stack_start: Define the depth with stacking beginning

      • stack_end: Define the depth with stacking ending

      • stack_val: The interval between adjacent nodes in depth

    • bin section: The shape of bin.

      • shape: Define the shape of bins. Only ‘rect’ or ‘circle’ is available

      • domperiod: the period in second of S wave used in calculating radius of a fresnel zone

      • witdth: The width of the profile.

      • bin_radius: Radius of each bin. If it is not defined, radius would be assumed following fresnel zone.

      • slide_val: Interval between adjacent bins

    • line section: Line location