Calculation of Theoretical Travel Time and Automatic Picking of Actual Travel Time in Seismic Data

: We used the ray tracing technique based on the IASP91 Earth model to calculate the travel times in order to identify the phases. This technique can calculate the travel times for the seismic phases in the conventional travel time tables. The waveform data received from the stations in the Guangxi area are selected for analysis and discussion. The outcomes of the numerical modeling and its use demonstrate that there is good agreement in terms of the absolute differences between the calculated and theoretical travel times from the ISAP91 tables. The relative residuals are determined directly from the actual arrival times picking during the correlation analysis, and the validity of the travel time method for picking seismic phases by correlation analysis can be demonstrated.


Introduction
For a determined model of the Earth's interior, it is possible to calculate from seismological methods what phases will occur at different epicentral distances and source depths according to the law of seismic wave propagation, and to determine the time they take to propagate from the source to the station (i.e., travel time), and, for an earthquake at a certain depth, the epicentral distance can be found from the arrival time difference between the subsequent phase and the first arrival phase [1].Seismic phase identification is a very important task in seismology.Seismic phase analysis is a powerful tool for studying seismic activity, source physics, earthquake engineering, and the internal physics of the Earth, and it is a fundamental tool for determining the three elements of an earthquake (time, epicenter and magnitude).
In 1935, Bullen and Jeffreys published the paper "Propagation times of seismic waves".In 1940, they collaborated on completing the Jeffreys-Bullen time table, referred to as the "J-B time table", which included most of the possible seismic waves on Earth [2].In the early 1980s, Dziewonski et al. obtained the Preliminary Reference Earth Model, a new generation of global one-dimensional seismic time model, which is the most accurate model compared to the J-B model.PREM is widely used in seismic research, especially in earthquake localization [3].Later, in order to meet the needs of global earthquake studies, the International Association of Seismology and Geophysics of the Earth's Interior (IASPEI), based on the preliminary Earth reference model, established an Earth model reflecting the global average properties in 1991, called the IASP91 model, and based on this model [4,5].Numerous ray tracing methods have been updated and improved as a result of the ongoing research and study of ray tracing techniques by various experts [6][7][8][9][10][11][12].
With the increasing number of seismic stations around the world, the seismic monitoring capability is increasing, but the consequent data generation is increasing exponentially, and how to process the massive seismic observation data quickly and effectively is a key problem faced.Seismic tomography has been widely applied to study the velocity structure of the Earth's interior and has achieved fruitful geological results [13][14][15][16].At present, the data used for laminar imaging of remote seismic body waves are usually the relative travel time residuals of the first arrival phase [15].The methods for obtaining first arrival times from waveform data mainly include manual identification and automatic computer detection methods.Although the manual identification method is simple to implement, it cannot meet the accuracy requirements for seismic waveform data with low signal-to-noise ratio, and it has a limited application because of its low efficiency and labor intensity; whereas the automatic computer detection method has become more and more popular among researchers due to its high accuracy and efficiency, such as the Akaike Information Criterion (AIC) arrival detection technique [17,18], the relative and absolute travel time difference simulated annealing method [19] and multiscale wavelet transform method [20].The short term average/long term average (STA/LTA) technique proposed by Allen R. V. (1978) [21] was used to perform automated picking.da Silva and Corso proposed a robust method based on instantaneous-spectral Shannon entropy for capturing microseismic events in noisy environments without data preprocessing [22].Wu et al. proposed an improved multi-channel cross-correlation method [23].Work on automatic picking through machine learning based technology is reported in [24,25].Russell B. proposed a research using machine learning and geophysical inversion.The author proposed that machine learning can be used for large amounts of data as an alternative to geophysical inversion [26].
Molyneux and Schmitt proposed the direct correlation method, but the study of the direct correlation method was carried out in a high time resolution, low noise signal obtained in the laboratory, so the scope of its use needs further study [27].VanDecar and Crosson proposed a multi-channel cross-correlation technique (MCCC) to achieve automatic acquisition of first arrival data, which not only greatly reduces human capital, but also improves data accuracy [28].Therefore, the multi-channel cross-correlation technique has been widely used in the field of seismic time-dependent data processing.However, the biggest drawback of the MCCC method is that the system of overdetermined equations must be solved when calculating the absolute travel time data, which inevitably reduces the accuracy of the travel time data.Then, the best absolute and relative travel time data are obtained by iteratively using a selected phase arrival time as the reference time [29].
In this study, the theoretical travel times of different seismic phases propagating in the Earth's interior are numerically simulated using a ray-tracing technique based on the IASP91 Earth velocity model.The measured travel times are quickly picked using a modified cross-correlation method, and the relative travel time residuals can be determined directly.The whole calculation process is completed automatically by the computer program, which not only saves human resources, but also improves the working efficiency.The waveform data received by the station in Guangxi area is selected as an example for analysis and discussion.

Seismic Phase Theory Travel Time Calculation
The IASP91 Earth model is shown in Figure 1 and Table 1.Note: x = r/R, where r is the radius from the Earth center, R is the Earth's radius (R = 6371.00km).
Using the IASP91 standard Earth velocity model for ray tracing [30], seismic waves propagate in a three-dimensional medium satisfying the equation function equation τ is the travel time of the seismic wave propagating from the source S to any point (x, y, z) in space.v is the velocity of the medium at any point in space.
Equation (1) can be expressed in the Cartesian coordinate system as Now, a variable s is introduced, which denotes the arc length of the ray as the elastic wave propagates.Since the ray is orthogonal to the wave front t = τ (x, y, z), the derivative of the ray at the wave front thus has Any spatial position (x, y, z) during the propagation of the elastic wave can be described by a parametric equation with the ray length s as a variable, i.e., x = x(s), y = y(s), z = z(s).For simplicity, the spatial position of the ray, i.e., the ray trajectory, can be expressed by a vector equation as Since the ray is orthogonal to the wavefront through which it passes, it follows that the ray path  must be parallel to ∇τ, i.e.,
From Equations ( 1) and ( 5), we can get Moreover, since ds 2 = dx 2 + dy 2 + dz 2 and thus we have α = v, it is possible to write the differential equation for the ray (5) in the following representation: From Equation (3), we get Finding the derivative with respect to s for both sides of Equation ( 7) and taking into account Equation ( 8), we have It can be seen that Equation ( 9) is a differential equation for the ray path that does not depend on the phase function τ.This equation can be expressed as three differential equations corresponding to the three coordinate components: If three directional cosines are induced at a point on the ray, i.e., from Equations ( 3) and ( 11) two equations can be obtained: Meanwhile, the other three first-order differential equations for the angle can be obtained from Equations ( 10) and (11): Note that in the above equation we have θx + θy = 90°, and using Equation (3), we get Equations ( 12) and ( 14) constitute the kinematic ray-tracing system in three-dimensional space.
For the IASP91 Earth velocity model, the medium velocity at any point inside the Earth is only related to the radial radius of that point, so the numerical simulation of the seismic wave travel time is reduced to ray tracing in a two-dimensional medium.The velocity of the medium at any point in the model is only related to the radial radius of that point, so the IASP91 Earth velocity model can be regarded as a two-dimensional circle structure model, and the ray paths of most of the seismic phases can be simulated by the semicircular circle structure model, except for several multiple reflection phases.
In this case, assuming that the velocity depends only on the x and z directions, and θx = α, θy = π/2, θz = π/2 − α, Equation ( 14) can be simplified as: Let the coordinates of the starting point S of the ray be (x0, z0) and the coordinates of any point P on the ray be (x, z), then Equation ( 4) can be discretized as: 0 is the off-source angle of the ray at the point (Figure 2), vx, vz are the velocity components of the seismic wave propagating to point P in the x and z coordinate axes, respectively, the off-source angle of the ray at point P (x, z), and Δτ is the time step of ray tracing.For the IASP91 Earth velocity model, it is more efficient to use the time-increment method for test shooting.The specific method is as follows: the source point S(x0, z0) is placed on the X-negative axis of the two-dimensional circle model; the ray departs from the source with the off-source angle α0, changes the size of the off-source angle α0 from 0° to 180°, and calculates the coordinates (x, z) and the off-source angle α of the ray arriving at point P after the elapsed time according to Equation ( 16).The corresponding ray path out of the surface point is regarded as the position of the detection point until the ray reaches the receiving position, completing a ray tracing process, and recording the corresponding travel time at different epicentral distances for different phases.
According to Snell's law there are: At each point of the ray, the ray parameter is p. v0 denotes the velocity of the ray at the epicenter S.
In this work, based on a simplified ray-tracing system in a two-dimensional medium (2D model is composed of concentric circles), ray tracing is performed using the timeincrement method [31].When a ray enters an interface (the partition between two layers of homogeneous medium) at an off-source angle, the off-source angle corresponding to the ray after refraction should be calculated according to Snell's law before ray tracing.When a ray propagates in a layer of homogeneous medium, the off-source angle at each point is calculated according to Equation ( 16).This method is simple and convenient to adjust the speed of each layer directly.

Automatic Picking of the Seismic Phase Actual Travel Time
In the same seismic event, there is a certain degree of similarity between stations, and the degree of similarity can be judged by the magnitude of the correlation coefficient, so the correlation coefficient is used as a measure of the degree of correlation between two random variables.
Let two columns of random signals x(t) and s(t) be accepted by N stations, and the correlation coefficient between x(t) = {x1, x2,…, xn} and s(t) = {s1, s2,…, sn} is The range of the correlation coefficient definition's value for Ri is 0 to 1.When Ri = 1, it indicates that the two series are identical; when Ri = 0, it indicates that the two series have no connection at all.The higher the value of Ri, the more similar the two waveforms are to one another.
The traditional correlation method is to determine the first arrival point by finding the delay time of the first arrival.However, because of the influence of the similarity between seismic traces, the delay time finding is prone to great deviation, and even the delay time cannot be obtained.When the waveform similarity between the reference channel and the other seismic channels is good, a good pickup effect can be achieved by directly using the reference channel to find the correlation coefficient with each seismic channel.However, the waveform information received by the seismic recording instruments is affected by various factors such as different propagation paths, noise, and receiving instruments, and the similarity of the actual data waveforms obtained is not ideal, so the pickup effect of the first arrivals is less accurate.
Many scholars pick by the similarity between the reference tract and each tract, and the correlation coefficient is found by the reference tract and each tract, and the correlation coefficient discrete form The discrete signal is xij (i = 1, 2, 3,…, N, j = 1, 2, 3,…, Nx); sk (k = 1, 2, 3, …, Ns), t = (k − 1)Δt, k = 1, 2, 3, …, Nx, Δt is the seismic waveform's sampling interval.The arrival time of the specified seismic phase received by the station is given by t = (k − 1)Δt if the correlation coefficient between the seismic record of the ith station after moment t and the specified wavelet is maximized.The seismic wave travel time ti of the phase can then be determined using the seismic event's onset moment t0.In Figure 3, the time signal sequences x(t) and s(t) were created by artificially shifting the same signal.x(t) was acquired after s(t) was artificially shifted to the right for 10 s, at the period when the value of Ri was at its highest.The data used in near-sound tomography are absolute travel times, while in farsound tomography, the data used are relative travel time residuals.
Suppose that the ith earthquake is observed by N stations after its occurrence.Where the actual arrival time of the first arrival phase manually picking from the waveform recorded at the jth station is denoted as  , and its theoretical arrival time is denoted as  , then the travel time residual of this phase tij is The average walk time residual  is: The relative travel time residual rij is: Equation ( 19) can be used to quickly calculate the arrival time difference data of all N waveforms with respect to the jth waveform, a total of N (where the arrival time difference of the jth waveform with respect to itself is 0 s), denoted as  , ⋯  , ( ) ⋯ 0 ⋯  , ( ) ⋯  , .The theoretical arrival time of each waveform can be calculated from the theoretical arrival time of each waveform to the first phase, and the theoretical arrival time of the jth waveform can be calculated as  , ⋯  , ( ) ⋯ 0 ⋯  , ( ) ⋯  , , which Make Then That is, the relative travel time residuals of Equation (22).The whole process is completely automated by the computer program, which largely improves the efficiency and accuracy of data processing. in Equation (20) represents the actual arrival time picking manually, and  in Equation ( 23) represents the actual arrival time variable, which is not involved in the actual operation.

Application Examples
The distribution of seismic stations is mostly located in Guangxi region, and Tables 2 and 3 shows the seismic records of two earthquakes observed by 16 stations of the broadband mobile seismic network in Guangxi region (Figure 4) as an example (rotated), #1.The Alo Islands experienced a magnitude 6.5 earthquake on 4 November 2015, at 03:44:15 (UTC = 0), at a source depth of 20 km, far from Guangxi.The P phase should be the initial arrival wave phase of this earthquake at a distance greater than 3000 km.With a source depth of 599.35 km and a distance of more than 17,000 km from Guangxi, the magnitude 6.7 earthquake in Brazil #2 occurred at 05:45:18 (UTC = 0) on 26 November 2015.The PKIKP phase should be the first arrival wave phase of this earthquake.The theoretical arrival time of the seismic phase was calculated according to the IASP91 Earth velocity model parameter table (Figures 5 and 6).The selected reference wavelet (SB11 station is used as the reference wavelet in Figure 7; SD21 is used as the reference wavelet in Figure 8).Use a shorter wavelet wherever possible to ensure the accuracy of the seismic phase picking.In general, the wavelet should not be longer than 15 s to prevent the wavelet from incorporating another seismic phase and influencing the accuracy of the travel time picking.The window length of the selected wavelets of other stations is set according to the length of the reference wavelet, so that the wavelets of other stations can be filtered to the best wavelets in the window set by the parameters (which will only be correlated within the window of the control parameters) even in the environment of low signal-to-noise ratio.shows the correlation coefficient, the coordinate axis is the upper red coordinate, the black dashed line is the relative travel time residual, and the coordinate axis is the lower black coordinate in seconds; the station SD21 is selected as the reference wavelet, and its correlation coefficient is 1).
The middle plots of Figures 9 and 10 are actually the results after correcting the Pwave seismic phase and PKIKP seismic phase when actual at each station with relative travel time residuals.It can be shown that the usage of cross-correlation travel time picking works effectively when combined with the correlation coefficient curves.The imperfect alignment of the P-phase and PKIKP-phase in Figures 9 and 10 reflects the presence of lateral structures below the stations.In the middle panel, most of the stations are perfectly aligned, in the physical sense by subtracting the relative residuals, and in the mathematical sense by pulling the wavelets of all 16 stations to the line where they were actual.In the right panel, it can be seen that the wavelets of each station are well correlated.

Conclusions
Teleseismic tomography is playing an increasingly important role in studying the velocity structure of the deep Earth.The accuracy of the relative travel time residuals used directly determines the credibility of the imaging results.Therefore, it has become an urgent problem for seismologists to obtain high-precision relative travel time residuals automatically from the huge amount of teleseismic waveform data, especially those with poor signal-to-noise ratio.
In this study, based on the IASP91 Earth model, we used a shooting method with a non-meshing scheme to trace raypaths and calculate the travel times of global seismic phases.The orthogonal travel times predicted by simulation for each phase can closely match the Kennett calculated seismic wave travel times for the corresponding phases.A comparison with SAC has also been done and the error time is within the range of 0.001 s.
This technique can show the kinematic properties of seismic phases moving within the Earth.The method does not require gridding of the Earth model and is simple to operate, fast to calculate, and highly accurate.
The advantage of correlation analysis is its simplicity and speed.There is a corresponding parameter window limit in calculating the number of interrelationships between two waveforms, and only the best wavelets are screened within the window range to improve the recognition of effective seismic phases.In the application of the teleseismic

Figure 2 .
Figure 2. Schematic diagram of speed difference interface point recalculation.

Figure 3 .
Figure 3.The cross-correlation functions of two time series.

Figure 5 .
Figure 5. Theoretical arrival times of seismic records and P and S phases in the Alo Islands.

Figures 9
Figures 9 and 10 show the cross-correlation effects, and the correlation coefficients are used to further assess the accuracy of the seismic phase travel time picking.

Figure 7 .
Figure 7. Correlation results of earthquake records in the Alo Islands.

Figure 8 .
Figure 8. Cross-correlation results of Brazilian earthquake records.

Figure 9 .
Figure 9. Cross-correlation evaluation system of seismic records in the Alo Islands.(The (left) figure shows the waveform with the theoretical arrival time of 0 line; the (middle) figure shows the waveform with the measured arrival time of 0 line; the red solid line in the (right) figure shows the correlation coefficient.The coordinate axis is the upper red coordinate, the black dashed line is the relative travel time residual, and the coordinate axis is the lower black coordinate in seconds; the station SB11 is selected as the reference wavelet, and its correlation coefficient is 1).

Figure 10 .
Figure 10.Cross-correlation evaluation system of seismic records of Brazilian earthquakes.(The (left) figure shows the waveform with the theoretical arrival time of 0 line; the (middle) figure shows the waveform with the measured arrival time of 0 line; the red solid line in the (right) figureshows the correlation coefficient, the coordinate axis is the upper red coordinate, the black dashed line is the relative travel time residual, and the coordinate axis is the lower black coordinate in seconds; the station SD21 is selected as the reference wavelet, and its correlation coefficient is 1).

Table 1 .
The parameters of IASP91 Earth model.

Table 2 .
Parameters of seismic events (parameters from USGS).