SplitRFLab: A MATLAB GUI toolbox for receiver function analysis based on SplitLab

We add new modules for receiver function (RF) analysis in SplitLab toolbox, which includes the manual RF analysis module, automatic RF analysis and related quality control modules, and H - k stacking module. The updated toolbox (named SplitRFLab toolbox), especially its automatic RF analysis module, could calculate the RFs quickly and efﬁciently, which is very useful in RF analysis with huge amount of seismic data. China is now conducting the ChinArray project that plans to deploy thousands of portable stations across Chinese mainland. Our SplitRFLab toolbox may obtain reliable RF results quickly at the ﬁrst time, which provide essentially new constraint to the crustal and mantle structures.


Introduction
Seismology is becoming the most popular tool in studying the interior structure of the Earth thanks to the development of digital seismology and the rapidly accumulated seismic data in recent years. Receiver function (RF) analysis, proposed by Langston (1979), is one of routine tools to determine the one-dimensional structure of the crust and mantle. P-wave RF deconvolves the radial component with the vertical component and obtains the Ps phase converted at velocity discontinuities (e.g., Moho) in the Earth. The time delay between the Ps converted phase and direct P-wave is sensitive to the depth of the velocity discontinuities and the average S-wave velocity above it (Zhu and Kanamori 2000). The RF method is stable and efficient in determining the crustal and upper mantle structures (e.g., Dueker and Sheehan 1997;Kind et al. 2002), which makes it one of the most popular methods in seismology.
The routines of RF analysis include the following procedures: (a) selecting the proper events (e.g., epicentral distance, magnitude) from the earthquake catalog and matching them with seismic data in the local directory; (b) preprocessing the waveforms, such as removing average shift of amplitude, removing the linear trend, and filtering the waveforms with certain bandpass filters; (c) picking P-wave arrivals; and (d) calculating P-wave RFs by deconvolving the radial component with the vertical component. Zhu and Kanamori (2000) further proposed the H-k stacking method with a grid-search algorithm that determines the v P /v S ratio and Moho depth (H) with the obtained P-wave RFs at a station, which is very important in determining one-dimensional structure in the Earth.
Many open-source codes have been released for RF analysis in recent years, such as the CPS 330 system, the RF module under Python, and the process RFmatlab module under MATLAB. These codes generally import the preprocessed seismic data to calculate RFs from teleseismic radial and vertical components. However, users have to write many small codes to obtain the preprocessed data first, e.g., selecting teleseismic events, removing average and linear trend from the waveforms, filtering, and selecting phases in proper time windows. Thanks to the rapid development of computer science, more GUI-based seismological softwares bring much convenience and power efficiency in seismic data analysis, such as the AIMBAT for picking arrival times (Lou et al. 2013), the FuncLab toolbox for RF analysis (Eagar and Fouch 2012), and SplitLab toolbox for shearwave splitting analysis (Wüstefeld et al. 2008).
SplitLab toolbox (Wüstefeld et al. 2008) developed under MATLAB environment by French scientist is an open-source GUI toolbox for shear-wave splitting analysis. It first links SAC data in the local directory to the earthquake catalog automatically, and then cuts the waveform in proper time windows by calculating the theoretical arrival times with Taup Toolkit (Crotwell et al. 1999). Users may preprocess the waveforms (e.g., removing average and linear trends, filtering with certain filters) quickly with defined keyboard shortcuts, which are then used for shear-wave splitting analysis and other related purposes. The GUI module of the toolbox provides very convenient and fast way to estimate the results interactively and efficiently, which makes it one of the most popular tools in shear-wave splitting analysis. FuncLab (Eagar and Fouch 2012) is a similar GUI-based toolbox for RF analysis. However, it lacks some nice functions such as automatically matching the local data with the earthquake catalog and interactively checking the results. Therefore, in this study, we integrate the RF analysis modules in the SplitLab toolbox, which inherits many nice functions for quick data preprocessing and convenient estimation of the final results. We also add the H-k stacking module (Zhu and Kanamori 2000) that obtains the crustal thickness and the average v P /v S in the crust under a station.

SplitRFLab toolbox
SplitRFLab, as an improved version of SplitLab toolbox, is a GUI-based software for RF and some other seismic data analysis under MATLAB environment. The toolbox can automatically import local SAC data, calculate their RFs, and selects high-quality results with convenient keyboard shortcuts. Besides, it includes modules to plot the results and conduct H-k stacking directly.

Data import module
SplitLab (Wüstefeld et al. 2008) provides very nice GUI modules for uses to link local SAC data to the proper event lists (constrained with dates, magnitudes, epicentral distances, azimuths, focal depths, etc.). It then calculates the theoretical arrivals of different phases with Taup Toolkit, and cut the waveforms in the corresponding time windows for further shear-wave splitting or RF analysis.
SplitLab also provides windows to view the seismograms conveniently. Users can select the waveforms in certain time windows by direct clicking and dragging with mouse in the window, which are prepared for subsequent processing such as time-frequency analysis and plotting particle motions. Because these excellent modules are also necessary in RF analysis, we keep them in our new SplitRFLab toolbox and add new functions and modules accordingly.

Parameter setting module
SplitRFLab adds many new modules to original SplitLab toolbox for RF analysis (Fig. 1). These modules include RF parameter setting module, manual RF calculation module, automatic RF calculation and related quality control modules, and H-k stacking module. In the parameter setting module (Fig. 2), users should set the output directory of the results as well as many necessary parameters used in RF analysis [e.g., length of waveforms, Gauss parameter, cutoff value of signal-to-noise ratio (SNR)]. Once these parameters are set, users many calculate the RFs and plot the results by simple clicks.

Manual RF analysis module
SplitRFLab toolbox adds functions of RF analysis in the original seismogram-viewer window (Fig. 3a). Similar to the shear-wave splitting analysis, users select the waveforms in the time windows containing P-wave which are converted to the radial, transverse, and vertical components automatically. Then users apply either iterative (Ligorría and Ammon 1999) or water-level (Ammon 1991) deconvolution algorithms by simple clicks or keyboard shortcuts to obtain the corresponding RFs. The codes for RF analysis are modified from process RFmatlab toolbox (https:// github.com/iwbailey/processRFmatlab). Figure 3b shows examples of radial and transverse RFs obtained by the manual RF analysis module; users may choose to save the results to the defined directory.
There are many functions in the seismogram-viewer window of SplitRFLab (Fig. 3), including zooming the waveforms, filtering, time-frequency analysis, coordinate transformation, etc. Users may call these functions by both mouse click and defined keyboard shortcuts (Wüstefeld et al. 2008). Here, we add new shortcuts for the RF calculation in the seismogram-viewer window; when the waveforms are selected, pressing keyboard ''R'' button will call functions to calculate the RFs.

Automatic RF analysis and quality control module
Compared with the manual RF analysis that obtains an RF with some clicks in the seismogram-viewer window, it is more efficient to calculate the RF automatically. Because the RF analysis is not sensitive to the accurate P-wave arrivals, we may first calculate the theoretical arrivals with Taup Toolkit (Crotwell et al. 1999) and then select the waveforms in an expanded time window that include the P-wave phase. The RFs are therefore calculated automatically with the selected preprocessed waveforms. Since the noise contained in different waveforms has different frequency bands, we calculate the RFs of waveforms with 3rd order Butterworth filter in three different frequency bands and then select one optimal result among them. These bands are sufficient for most studies in P-wave receiver function analysis. However, users may change the desired frequency bands in ''calRFbat.m'' according to the quality of their data. We also calculate the SNR value of P-wave by comparing the energy of waveforms (10 s long) after and before the theoretical arrivals. We first set cut-off values of SNR for the radial and vertical components in the parameter setting module. Only the waveforms with SNR values larger than the cut-off values are processed further to calculate the RFs, which enhances the efficiency significantly. SplitRFLab adds a new window to plot the obtained RFs (Fig. 4), which are used to check and select the optimal result from the results of three different frequency bands. We defined keyboard shortcuts (Table 1) for users to select the optimal results quickly. For example, pressing keyboard ''q'' button in the window will save the result of the waveforms filtered in the frequency bands of 0.03-2.0 Hz. Users may also press keyboard ''s'' button to import and view the original waveforms quickly, which provides important constraints on the reliability of the obtained RFs. More keyboard shortcuts can be defined in the file ''rfKeyPress.m'' when necessary.

H-k stacking module
After calculating and saving the RFs of all events at a station, users may plot the RF section and save it in local directory by simple clicking at defined buttons. The SplitRFLab also provides window for H-k stacking; users set proper parameters and then press the button to calculate the v P /v S value and Moho depth (H) with a grid-search algorithm (Zhu and Kanamori 2000). We estimate the uncertainties of (v P /v S , H) pair following Eaton et al. (2006) and plot an ellipse in the figure to show them visually. All the results, including the figures and (v P /v S , H) values, are finally saved in the predefined directory.

Quick guide to SplitRFLab
The complete modules and functions of SplitRFLab toolbox are shown in Fig. 1, in which the red squares and red words in blue squares are newly added for RF analysis. Here, we provide an example of calculating the RFs and obtaining the crustal structure under a permanent station NJ2 (Fig. 5) with SplitRFLab toolbox.
We first link the local SAC data to the earthquake catalog under the instruction of SplitLab (Wüstefeld et al. 2008). Then in the panel ''Process'' we set the path to the directories where the cut waveforms, calculated RFs, and RF section are saved (Fig. 2). Clicking the button ''RF parameters'' will call a pop-up window where we set the values of Gauss parameter, method, and cut-off values of SNR that are necessary in RF calculation. Now click the button ''Do PRFs'', and the SplitRFLab toolbox will calculate the RFs of all events at the station NJ2. A pop-up bar appears and shows the process of the calculation. Those waveforms with SNR values smaller than the predefined cut-off values are generally ignored in the step, which accelerates the calculation. After the automatic RF analysis, we click the button ''View RF'' in panel ''View  Table 1 Database'', check every RF, and save good results with keyboard shortcuts (Table 1). At last, we click the button ''Plot RT'' and a figure showing the RF section plotted with a function of back azimuth will appear (Fig. 6).
The H-k stacking module makes it possible to obtain the (v P /v S , H) values quickly. We first set ranges of the parameters (e.g., v P /v S , H) in the grid-search algorithm as well as average P-wave velocity in the crust, weights of different phases in the stacking, and the path to the directory for the saved plots and results. Then click button ''Go stacking'', the SplitRFLab will call the H-k stacking module to obtain the energy map of different phases (i.e., Ps, PpPs, PsPs, & PpSs) (Fig. 7a-c) and their weight summary (Fig. 7d). The optimal (v P /v S , H) pair is also calculated automatically. We may choose to save the plot and result or not. Moreover, clicking the button ''Plot R'' will plot the radial RFs as a function of ray parameter (Fig. 7e). The theoretical arrivals of different phases with the optimal (v P /v S , H) model will be marked in the RF section (Fig. 7e).

Summary
We developed the SplitRFLab code, which is a MATLAB toolbox with friendly GUI, and added new RF analysis 30 60 90 NJ2 Fig. 5 The events used in RF analysis at station NJ2  Blue and white colors show strong and weak energy, respectively. The optimal (v P /v S , H) pair is marked as red lines while the ellipse denotes their uncertainties (Eaton et al. 2006). e Radial RFs at station NJ2 plotted as a function of ray parameters; the stacking RF is shown at the top. The three-dashed vertical lines denote the theoretical arrival times of the Ps, PpPs, and PpSs ? PsPs phases Earthq Sci (2016) 29(1):17-26 25