Logo

Stack PRFs by Common Conversion points (CCP) method

This note will introduce how to stack PRFs by 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

Befor the CCP stacking, we should prepare:

1.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.

2. 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.

3. A parameter file

The parameter file include parameters used in CCP stacking. the format should follow configparser as the following example:

[FileIO]
# Path to stations with RF sac files
rfpath = /path/to/RFresult

# Path to station list
stalist = /path/to/sta_all.lst

# Path to the lib of Ps ray parameters. 
# If it's empty the ray parameters of Ps would be assumed as that of P arrival
rayp_lib =

# Output data struct after time to depth
depthdat =  /Users/xumj/Researches/NETibet/ccp_result/RFdepth_3D.mat

# Output data struct after CCP stacking
stackfile = /Users/xumj/Researches/NETibet/Ordos_Process/stack_L5

# Station list used to stack
stack_sta_list = /Users/xumj/Researches/NETibet/Ordos_Process/sta_L5.lst

# Path to 1D velocity model
# If it's empty, useing IASP91 model
velmod =

[bin]
# The shape of bins
shape = rect

# period of S wave (to assuming the radius of fresnel zone)
domperiod = 5

# Width of the profile in km
width = 100

# Radius of bins in km
# If it's empty, radius would be assumed following fresnel zone
bin_radius =

# sliding interval of bins in km 
slid_val = 5

[line]
# Coordinate of two end points for the profile
profile_lat1 = 40.82
profile_lon1 = 103.451
profile_lat2 = 36.33
profile_lon2 = 108.817

[depth]
# Max depth you wanna convert to
dep_end = 800

# depth interval in converting time to depth
dep_val = 1

[stack]
# Stack RFs from <stack_start> km to <stack_end> km with interval of <stack_val> km
stack_start = 0
stack_end = 150
stack_val = 1

4. (Optional) A lib of Ps ray parameters

If you would stack PRFs in a great 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 DIS_STR -e DEP_STR [-l LAY_STR] [-o OUT_PATH]

Gen a ray parameter lib for Pds phases

optional arguments:
  -h, --help   show this help message and exit
  -d DIS_STR   Distance range as 'min_dis/max_dis/interval'
  -e DEP_STR   Event depth range as 'min_dep/max_dep/interval'
  -l LAY_STR   layers range as 'min_layer/man_layer'
  -o OUT_PATH  Out path to Pds ray parameter lib

5. (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 VEL3DPATH] 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 VEL3DPATH  Path to 3d vel model in npz file

The output struct 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] [-l LINE_LOCA] [-s STALIST] [-o OUTPATH] [-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
  -l LINE_LOCA  Location of the profile <lat1>/<lon1>/<lat2>/<lon2>
  -s STALIST    Path to station list
  -o OUTPATH    Path to output data
  -t            Output as text file

Required arguments

  • cfg_file: the parameter file introduced above. Meanwhile, some options are required:
    • FileIO: Allow to set 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 witdth of the profile. the stations that meets the condition would be project to this proflie and write to this stalist.
    • stack section: Define the nodes in depth
      • stack_start: Define the depth with stacking begining
      • stack_end: Defane 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 avalible
      • 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.
      • slid_val: Interval between adjacent bins line: Line location