Next Article in Journal
A Comparative Assessment of JVM Frameworks to Develop Microservices
Next Article in Special Issue
Estimating Liquefaction Susceptibility Using Machine Learning Algorithms with a Case of Metro Manila, Philippines
Previous Article in Journal
Explore Long-Range Context Features for Speaker Verification
Previous Article in Special Issue
An Interactive System Based on the IASP91 Earth Model for Earthquake Data Processing
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

College of Earth Sciences, Guilin University of Technology, Guilin 541004, China
*
Authors to whom correspondence should be addressed.
Appl. Sci. 2023, 13(3), 1341; https://doi.org/10.3390/app13031341
Submission received: 6 December 2022 / Revised: 12 January 2023 / Accepted: 16 January 2023 / Published: 19 January 2023
(This article belongs to the Special Issue Big Data in Seismology: Methods and Applications)

Abstract

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

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

2. Seismic Phase Theory Travel Time Calculation

The IASP91 Earth model is shown in Figure 1 and Table 1.
Using the IASP91 standard Earth velocity model for ray tracing [30], seismic waves propagate in a three-dimensional medium satisfying the equation function equation
τ 2 = 1 v 2  
τ 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
τ x 2 + τ y 2 + τ z 2 = 1 v 2  
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
d τ d s = 1 v
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
x ˜ = x ˜ s
Since the ray is orthogonal to the wavefront through which it passes, it follows that the ray path x ˜ must be parallel to ∇τ, i.e.,
d x ˜ d s = α τ
where α is the scaling factor.
From Equations (1) and (5), we can get
τ 2 = 1 α 2 d x ˜ d s 2 = 1 α 2 d x d s 2 + d y d s 2 + d z d s 2 = 1 v 2  
Moreover, since ds2 = dx2 + dy2 + dz2 and thus we have α = v, it is possible to write the differential equation for the ray (5) in the following representation:
1 v d x ˜ d s = τ
From Equation (3), we get
d τ d s = 1 v
Finding the derivative with respect to s for both sides of Equation (7) and taking into account Equation (8), we have
d d s 1 v d x ˜ d s = d d s τ = 1 v
i.e.,
d d s 1 v d x ˜ d s = 1 v
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:
d d s 1 v d x d s = x 1 v d d s 1 v d y d s = y 1 v d d s 1 v d z d s = z 1 v
If three directional cosines are induced at a point on the ray, i.e.,
c o s θ x = x s c o s θ y = y s c o s θ z = z s
from Equations (3) and (11) two equations can be obtained:
x τ = v c o s θ x y τ = v c o s θ y z τ = v c o s θ z
Meanwhile, the other three first-order differential equations for the angle can be obtained from Equations (10) and (11):
d d s 1 v c o s θ x = x 1 v d d s 1 v c o s θ y = y 1 v d d s 1 v c o s θ z = z 1 v   i . e . ,   d d s 1 v c o s θ x 1 v s i n θ x d θ x d s = x 1 v d d s 1 v c o s θ y 1 v s i n θ y d θ y d s = y 1 v d d s 1 v c o s θ z 1 v s i n θ z d θ z d s = z 1 v
Note that in the above equation we have θx + θy = 90°, and using Equation (3), we get
θ x τ = v x s i n θ x v y c o t θ x c o s θ y v z c o t θ x c o s θ z θ y τ = v y s i n θ y v x c o t θ x c o s θ y v z c o t θ y c o s θ z θ z τ = v z s i n θ z v x c o t θ x c o s θ z v y c o t θ x c o s θ y
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:
x τ = v c o s α   z τ = v s i n α   α τ = v x s i n α v z c o s α
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:
x = x 0 + v s i n α 0 · Δ τ   z = z 0 + v c o s α 0 · Δ τ   α = α 0 v x c o s α 0 v z s i n α 0 · Δ τ
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:
s i n α 0 v 0 = s i n α v x , z = p  
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 time-increment 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.

3. 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
R i τ = x i t + τ s t x i 2 t s 2 t
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
R i k = j = k k + N s 1 x i j s j j = k k + N s 1 x i j 2 j = k k + N s s j 2
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 far-sound 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 T i j o b s , and its theoretical arrival time is denoted as T i j c a l , then the travel time residual of this phase tij is
t i j = T i j o b s T i j c a l
The average walk time residual   t i ¯ is:
t i ¯ = 1 N j = 1 N t i j
The relative travel time residual rij is:
r i j = t i j t i ¯
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 Y i , j 1 o b s     Y i , j j 1 o b s   0     Y i , j j + 1 o b s     Y i , j N o b s . 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 Y i , j 1 cal     Y i , j j 1 c a l   0     Y i , j j + 1 c a l     Y i , j N c a l , which
  Y i , j k o b s = T i j o b s T i k o b s
Y i , j k c a l = T i j c a l T i k c a l
Make
  Y i , j k o b s ¯ = 1 N k = 1 N Y i , j k o b s
Y i , j k c a l ¯ = 1 N k = 1 N Y i , j k o b s
Then
y i j = Y i , j k o b s ¯ Y i , j k c a l ¯ = 1 N k = 1 N Y i , j k o b s Y i , j k c a l = 1 N k = 1 N Y i , j o b s Y i , j c a l Y i , k o b s Y i , k c a l = 1 N k = 1 N t i j t i k = t i j 1 N k = 1 N t i k = t i j t i ¯ = r i j
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. T i j o b s in Equation (20) represents the actual arrival time picking manually, and T i j o b s in Equation (23) represents the actual arrival time variable, which is not involved in the actual operation.

4. Application Examples

The distribution of seismic stations is mostly located in Guangxi region, and Table 2 and Table 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 (Figure 5 and Figure 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.
Figure 9 and Figure 10 show the cross-correlation effects, and the correlation coefficients are used to further assess the accuracy of the seismic phase travel time picking.
The middle plots of Figure 9 and Figure 10 are actually the results after correcting the P-wave 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 Figure 9 and Figure 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.

5. 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 cross-correlation evaluation system, the theoretical arrival times and actual arrival times calculated using the IASP91 Earth velocity model can roughly align with the initial waveforms, and the relative residuals are determined directly from the actual arrival times picked up during correlation analysis without any inversion or reference time information. The effectiveness of the travel time method for picking the seismic phase by correlation analysis can be demonstrated.

Author Contributions

Conceptualization, W.G. and S.P.; methodology, W.G.; software, S.Y.; validation, Y.W., S.P. and S.Y.; formal analysis, W.G.; writing—original draft preparation, W.G. and S.P.; writing—review and editing, L.L. and Y.W.; supervision, L.Y. and Y.Y. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by the National Natural Science Foundation of China (No. 42062015), Innovation Program for College Graduates of Guangxi (No. YCSW2022312), and the Natural Science Foundation of Guangxi (No. 2018AD19204).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data available in a publicly accessible repository (https://earthquake.usgs.gov/) (accessed on 1 December 2022).

Acknowledgments

Thanks to the USGS website for the data.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Huang, Z. Development of Automatic Earthquake Quick Report in China. Earthq. Res. China 2020, 34, 219–226. [Google Scholar]
  2. Jeffreys, H.; Bullen, K.E. Seismology Tables; British Association of Seismological Investigations: London, UK, 1940. [Google Scholar]
  3. Dziewonski, A.M.; Anderson, D.L. Preliminary reference Earth model. Phys. Earth Planet. Inter. 1981, 25, 297–356. [Google Scholar] [CrossRef]
  4. Kenntt, B.L.N.; Engdahl, E.R. Travel times for global earthquake location and phase association. Geophys. J. Int. 1991, 105, 429–465. [Google Scholar] [CrossRef] [Green Version]
  5. Engdahl, E.R.; Van der Hilst, R.; Buland, R. Global teleseismic earthquake relocation with improved travel times and procedures for depth determination. Bull. Seismol. Soc. Am. 1998, 88, 722–743. [Google Scholar] [CrossRef]
  6. Steck, L.K.; Thurber, C.H.; Fehler, M.C. Crust and upper mantle P wave velocity structure beneath Valles caldera, New Mexico: Results from the Jemez teleseismic tomography experiment. J. Geophys. Res. Solid Earth 1998, 103, 24301–24320. [Google Scholar] [CrossRef]
  7. Bijwaard, H.; Spakman, W. Fast kinematic ray tracing of first-and later-arriving global seismic phases. Geophys. J. Int. 1999, 139, 359–369. [Google Scholar] [CrossRef] [Green Version]
  8. Keyser, M.; Ritter, J.R.R.; Jordan, M. 3D shear-wave velocity structure of the Eifel plume, Germany. Earth Planet. Sci. Lett. 2002, 203, 59–82. [Google Scholar] [CrossRef]
  9. Zhao, D.; Lei, J. Seismic ray path variations in a 3D global velocity model. Phys. Earth Planet. Inter. 2004, 141, 153–166. [Google Scholar] [CrossRef]
  10. Wang, J.Y. Inverse Theory in Geophysics, 2nd ed.; Higher Education Press: Beijing, China, 2002. [Google Scholar]
  11. Bai, C.; Huang, G.; Zhao, R. 2-D/3-D irregular shortest-path ray tracing for multiple arrivals and its applications. Geophys. J. Int. 2010, 183, 1596–1612. [Google Scholar] [CrossRef] [Green Version]
  12. Huang, G.J.; Bai, C.Y.; Greenhalgh, S. Fast and accurate global multiphase arrival tracking: The irregular shortest-path method in a 3-D spherical earth model. Geophys. J. Int. 2013, 194, 1878–1892. [Google Scholar] [CrossRef] [Green Version]
  13. Guo, B.; Liu, Q.Y.; Chen, J.H.; Zhao, D.P.; Li, S.C.; Lai, Y.G. Seismic tomography of the Seismic tomography of the crust and upper mantle structure underneath the Chinese Tianshan. Chin. J. Geophys. 2006, 49, 1693–1700. (In Chinese) [Google Scholar] [CrossRef]
  14. Zhao, D.; Hasegawa, A.; Horiuchi, S. Tomographic imaging of P and S wave velocity structure beneath northeastern Japan. J. Geophys. Res. Solid Earth 1992, 97, 19909–19928. [Google Scholar] [CrossRef]
  15. Zhao, D.; Hasegawa, A.; Kanamori, H. Deep structure of Japan subduction zone as derived from local, regional, and teleseismic events. J. Geophys. Res. Solid Earth 1994, 99, 22313–22329. [Google Scholar] [CrossRef] [Green Version]
  16. Huang, J.; Zhao, D. High-resolution mantle tomography of China and surrounding regions. J. Geophys. Res. Solid Earth 2006, 111, B09305. [Google Scholar] [CrossRef]
  17. Sleeman, R.; Van Eck, T. Robust automatic P-phase picking: An on-line implementation in the analysis of broadband seismogram recordings. Phys. Earth Planet. Inter. 1999, 113, 265–275. [Google Scholar] [CrossRef]
  18. Wang, J.; Chen, J.H.; Liu, Q.Y.; Li, S.C.; Guo, B. Automatic onset phase picking for portable seismic array observation. Acta Seismol. Sin. 2006, 28, 42–51. (In Chinese) [Google Scholar] [CrossRef]
  19. Bear, L.K.; Pavlis, G.L. Multi-channel estimation of time residuals from broadband seismic data using multi-wavelets. Bull. Seismol. Soc. Am. 1999, 89, 681–692. [Google Scholar] [CrossRef]
  20. Chevrot, S. Optimal measurement of relative and absolute delay times by simulated annealing. Geophys. J. Int. 2002, 151, 164–171. [Google Scholar] [CrossRef] [Green Version]
  21. Allen, R.V. Automatic earthquake recognition and timing from single traces. Bull. Seismol. Soc. Am. 1978, 68, 1521–1532. [Google Scholar] [CrossRef]
  22. da Silva, S.L.E.F.; Corso, G. Microseismic event detection in noisy environments with instantaneous spectral Shannon entropy. Phys. Rev. E 2022, 106, 014133. [Google Scholar] [CrossRef]
  23. Wu, H.; Xiao, W.; Ren, H. Automatic Time Picking for Weak Seismic Phase in the Strong Noise and Interference Environment: An Hybrid Method Based on Array Similarity. Sensors 2022, 22, 9924. [Google Scholar] [CrossRef] [PubMed]
  24. Dimililer, K.; Dindar, H.; Al-Turjman, F. Deep learning, machine learning and internet of things in geophysical engineering applications: An overview. Microprocess. Microsyst. 2021, 80, 103613. [Google Scholar] [CrossRef]
  25. Wang, J.; Xiao, Z.; Liu, C.; Zhao, D.; Yao, Z. Deep learning for picking seismic arrival times. J. Geophys. Res. Solid Earth 2019, 124, 6612–6624. [Google Scholar] [CrossRef]
  26. Russell, B. Machine learning and geophysical inversion—A numerical study. Lead. Edge 2019, 38, 512–519. [Google Scholar] [CrossRef]
  27. Molyneux, J.B.; Schmitt, D.R. First-break timing: Arrival onset times by direct correlation. Geophysics 1999, 64, 1492–1501. [Google Scholar] [CrossRef] [Green Version]
  28. VanDecar, J.C.; Crosson, R.S. Determination of teleseismic relative phase arrival times using multi-channel cross-correlation and least squares. Bull. Seismol. Soc. Am. 1990, 80, 150–169. [Google Scholar]
  29. Rawlinson, N.; Kennett, B.L.N. Rapid estimation of relative and absolute delay times across a network by adaptive stacking. Geophys. J. Int. 2004, 157, 332–340. [Google Scholar] [CrossRef]
  30. Jiang, C.; Wang, Y.; Xiong, B.; Ren, Q.; Hu, J.; Gao, W.; Tian, Y.; Xi, Z. Numerical modeling of global seismic phases and its application in seismic phase identification. Earthq. Sci. 2019, 32, 72–79. [Google Scholar] [CrossRef]
  31. Boisse, P.H.; Bussy, P.; Ladeveze, P. A new approach in non-linear mechanics: The large time increment method. Int. J. Numer. Methods Eng. 1990, 29, 647–663. [Google Scholar] [CrossRef]
Figure 1. IASP91 Earth model.
Figure 1. IASP91 Earth model.
Applsci 13 01341 g001
Figure 2. Schematic diagram of speed difference interface point recalculation.
Figure 2. Schematic diagram of speed difference interface point recalculation.
Applsci 13 01341 g002
Figure 3. The cross-correlation functions of two time series.
Figure 3. The cross-correlation functions of two time series.
Applsci 13 01341 g003
Figure 4. Map of stations.
Figure 4. Map of stations.
Applsci 13 01341 g004
Figure 5. Theoretical arrival times of seismic records and P and S phases in the Alo Islands.
Figure 5. Theoretical arrival times of seismic records and P and S phases in the Alo Islands.
Applsci 13 01341 g005
Figure 6. Brazilian earthquake records and theoretical arrival times of PKIKP, PKPab, SKIKP, PKIKS, PP, SKIKS, and PPP seismic phases.
Figure 6. Brazilian earthquake records and theoretical arrival times of PKIKP, PKPab, SKIKP, PKIKS, PP, SKIKS, and PPP seismic phases.
Applsci 13 01341 g006
Figure 7. Correlation results of earthquake records in the Alo Islands.
Figure 7. Correlation results of earthquake records in the Alo Islands.
Applsci 13 01341 g007
Figure 8. Cross-correlation results of Brazilian earthquake records.
Figure 8. Cross-correlation results of Brazilian earthquake records.
Applsci 13 01341 g008
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 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).
Applsci 13 01341 g009
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) 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 SD21 is selected as the reference wavelet, and its correlation coefficient is 1).
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) 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 SD21 is selected as the reference wavelet, and its correlation coefficient is 1).
Applsci 13 01341 g010
Table 1. The parameters of IASP91 Earth model.
Table 1. The parameters of IASP91 Earth model.
Depth (km)Radius (km)VP (km/s)VS (km/s)
0.0~20.06351~63715.83.36
20.0~35.06336~63516.53.75
35.0~120.06251~63368.78541 − 0.749530x6.706231 − 2.248585x
120.0~210.06161~625125.41389 − 17.69722x5.750200 − 1.274200x
210.0~410.05961~616130.78765 − 23.25415x15.242130 − 11.085520x
410.0~660.05711~596129.38896 − 21.40656x17.707320 − 13.506520x
660.0~760.05611~571125.96984 − 16.93412x20.768900 − 16.531470x
760.0~2740.03631~561125.14860 − 41.15380x + 51.99320x2 − 26.60830x312.930300 − 21.259000x + 27.89880 x2 − 14.10800 x3
2740.0~2889.03482~363114.49470 − 1.47089x8.166160 − 1.582060x
2889.0~5153.91217.1~348210.03904 + 3.75665x − 13.67046 x20
5153.9~6371.00~1217.111.24094 − 4.09689 x23.564540 − 3.452410 x2
Note: x = r/R, where r is the radius from the Earth center, R is the Earth’s radius (R = 6371.00 km).
Table 2. Parameters of seismic events (parameters from USGS).
Table 2. Parameters of seismic events (parameters from USGS).
EventTime (UTC)LongitudeLatitudeDepth (km)MagnitudePhase
#14 November 2015 3:44:15124.88−8.3420.006.5P
#226 November 2015 5:45:18−71.29−9.19599.356.7PKIKP
Table 3. Relative travel time difference and correlation coefficient.
Table 3. Relative travel time difference and correlation coefficient.
StatLat.Lon.Alt. (m)#1#2
DistanceTobsTcalTdeltCoefDistanceTobsTcalTdeltCoef
YK2121.9943 111.0042 97 33.193 394.642 395.242 −0.60.933 166.948 1138.193 1137.163 1.03 0.949
FS2122.0011 109.9945 84 33.617 398.328 398.928 −0.60.909 167.069 1138.244 1137.254 0.99 0.945
YA1123.0113 111.9918 39 33.739 399.389 399.989 −0.60.990 165.767 1137.350 1136.240 1.11 0.943
GL2122.0069 109.4934 65 33.838 400.549 400.849 −0.30.941 167.100 1138.087 1137.277 0.81 0.957
NP2121.9902 108.9901 2 34.046 402.453 402.653 −0.20.976 167.136 1137.884 1137.304 0.58 0.943
ZC1123.0118 109.0052 86 34.937 410.847 410.347 0.50.986 166.114 1137.632 1136.522 1.11 0.918
ZP2124.0041 111.0215 120 35.007 411.052 410.952 0.10.971 164.962 1136.745 1135.575 1.17 0.964
SB1123.5074 109.4913 135 35.169 412.608 412.340 0.2691.000 165.602 1136.938 1136.108 0.83 0.918
BB2124.5018 111.4985 153 35.287 413.154 413.354 −0.20.941 164.399 1136.211 1135.081 1.13 0.959
JX2124.0076 109.9857 183 35.414 414.745 414.445 0.30.953 165.072 1136.518 1135.668 0.85 0.976
XP2124.4994 110.4897 194 35.660 416.958 416.558 0.40.987 164.536 1136.256 1135.206 1.05 0.900
JY2125.0271 111.0226 339 35.941 419.264 418.964 0.30.986 163.950 1135.772 1134.682 1.09 0.946
YF2125.0157 110.0182 138 36.309 422.412 422.112 0.30.976 164.065 1135.537 1134.787 0.75 0.979
SD2125.0039 109.5137 358 36.497 424.116 423.716 0.40.988 164.107 1135.508 1134.825 0.68 1.000
LC2125.0046 108.9980 137 36.707 425.899 425.499 0.40.990 164.122 1135.398 1134.839 0.56 0.949
SJ2125.9423 109.5843 240 37.314 431.058 430.658 0.40.986 163.166 1134.542 1133.952 0.59 0.968
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Gao, W.; Wang, Y.; Yang, Y.; Peng, S.; Yu, S.; Liu, L.; Yan, L. Calculation of Theoretical Travel Time and Automatic Picking of Actual Travel Time in Seismic Data. Appl. Sci. 2023, 13, 1341. https://doi.org/10.3390/app13031341

AMA Style

Gao W, Wang Y, Yang Y, Peng S, Yu S, Liu L, Yan L. Calculation of Theoretical Travel Time and Automatic Picking of Actual Travel Time in Seismic Data. Applied Sciences. 2023; 13(3):1341. https://doi.org/10.3390/app13031341

Chicago/Turabian Style

Gao, Wenqi, Youxue Wang, Yang Yang, Sanxi Peng, Songping Yu, Lu Liu, and Lei Yan. 2023. "Calculation of Theoretical Travel Time and Automatic Picking of Actual Travel Time in Seismic Data" Applied Sciences 13, no. 3: 1341. https://doi.org/10.3390/app13031341

APA Style

Gao, W., Wang, Y., Yang, Y., Peng, S., Yu, S., Liu, L., & Yan, L. (2023). Calculation of Theoretical Travel Time and Automatic Picking of Actual Travel Time in Seismic Data. Applied Sciences, 13(3), 1341. https://doi.org/10.3390/app13031341

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop