Seismic imaging travel time data processing method and device

Abstract

The present invention provides a method and apparatus for processing seismic imaging travel time data. The method comprises: establishing a seismic velocity structure model of a target work area based on seismic imaging data of the target work area; establishing an eikonal equation for the seismic imaging data; and determining the earthquake first arrival times of the seismic imaging travel time data based on the seismic velocity structure model, the eikonal equation, and pre-established upwind difference format data. The method and apparatus for processing seismic imaging travel time data provided by the present invention can replace manual processing of seismic imaging travel time data, achieving fully automated processing of seismic imaging travel time data and offering the advantages of high efficiency and excellent accuracy.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic data processing, and in particular to a method and device for processing seismic imaging travel time data. Background Art

[0002] In the existing technology, the use of geophysical fields to detect the Earth's internal structure is considered the most effective means of studying the Earth's internal structure. Specifically, the Earth's internal structure can be described through seismic imaging technology. It can be understood that seismic imaging uses the elastic wave field excited by earthquakes to image the Earth's internal structure. Seismic imaging can image parameters such as the velocity structure, anisotropy structure, attenuation structure, and interface structure of the underground medium. Today, popular seismic imaging methods can be mainly divided into two categories. One is adjoint imaging based on the wave equation; the other is ray travel-time imaging.

[0003] For the first type of method, adjoint imaging based on the wave equation performs waveform cross-correlation between the forward wavefield from the source and the reverse wavefield from the station. It then constructs the Frechét derivatives of the error functional with respect to the model parameters. Once these Frechét derivatives are obtained, numerical optimization methods such as the steepest descent method, the conjugate gradient method, Newton's method, and the BFGS method are used to iteratively optimize and gradually obtain the target model. Adjoint imaging based on the wave equation uses the wave method to construct the Frechét derivatives. Because this method incorporates both kinematic (e.g., travel time) and dynamic (e.g., amplitude and finite frequency) information about seismic waves, it offers high imaging resolution (maximum resolution can reach half the shortest wavelength). However, both the numerical solution of the source wavefield and the solution of the reverse wavefield from the station require numerical wave equations, which require significant computational resources. Therefore, adjoint imaging using the wave equation is costly and inefficient.

[0004] Regarding the second type of method, ray traveltime imaging offers high imaging efficiency. Under high-frequency approximations, seismic wave propagation can be considered as ray propagation. Ray traveltime imaging primarily constructs the Frechét derivatives of traveltime functionals with respect to model parameters by analyzing the sensitivity of ray traveltimes to the model. Calculating seismic ray paths is a key step in ray traveltime imaging. There are two main approaches for calculating ray paths: the pseudo-bending method and the eikonal equation solution method. The pseudo-bending method perturbs the ray path between the source and the station in a specific manner, continuously correcting the ray path to obtain the final calculated ray path. The pseudo-bending method involves calculating the spatial first-order derivative of the model velocity, which requires a smooth model. However, the real earth medium often has discontinuities, making the pseudo-bending method difficult to apply to seismic imaging problems in complex media. The eikonal equation solution method considers the impact of medium heterogeneity on traveltimes and is applicable to any heterogeneous medium model, demonstrating high model adaptability. Numerical solutions to the eikonal equation are computationally efficient and do not require the use of massively parallel computers. Eikonal equation imaging offers a good balance between imaging accuracy and cost.

[0005] When using the eikonal equation for imaging, the inverted data object is traveltime residuals, which can be obtained by calculating the difference between observed traveltimes and synthesized traveltimes. Synthesized traveltimes are relatively easy to obtain. After solving the eikonal equation to obtain traveltimes at spatial grid points, the synthesized traveltimes at the station can be obtained through spatial interpolation. However, obtaining observed traveltimes requires manual calibration of the actual data. The accuracy of traveltime calibration in the actual data directly affects imaging accuracy. Observed traveltimes are typically calibrated manually by trained personnel. While manual identification offers high accuracy and reliable imaging results, it is labor-intensive and slow. Currently, the number of earthquake traveltimes requiring calibration for seismic imaging reaches hundreds of thousands or even millions. Manual identification alone is no longer sufficient for processing the massive amounts of seismic data.

[0006] In summary, there is an urgent need for an automatic processing method for seismic travel time data to improve the data processing capabilities in the field of seismic imaging. Summary of the Invention

[0007] In response to the problems in the prior art, the seismic imaging travel time data processing method and device provided by the present invention can replace manual processing of seismic imaging travel time data, realize full automation of seismic imaging travel time data processing, and have the advantages of high efficiency and excellent precision.

[0008] In order to solve the above technical problems, the present invention provides the following technical solutions:

[0009] In a first aspect, the present invention provides a method for processing seismic imaging travel time data, comprising:

[0010] Establishing a seismic velocity structure model of the target work area based on seismic imaging data of the target work area;

[0011] Establishing an eikonal equation for the seismic imaging data;

[0012] The earthquake first arrival time of the seismic imaging travel time data is determined according to the seismic velocity structure model, the eikonal equation and pre-established upwind difference format data.

[0013] In one embodiment, establishing the eikonal equation for the seismic imaging data includes:

[0014] The eikonal equation is established based on the spatial travel time field, spatial gradient operator, medium slowness of the seismic imaging data and basic data of the target work area.

[0015] In one embodiment, the method for establishing the upwind difference format data includes:

[0016] Determining an upwind difference operator according to the spatial travel time field of the seismic imaging data and its grid characterization parameters;

[0017] The upwind difference format data is established according to the upwind difference operator.

[0018] In one embodiment, determining the earthquake first arrival time of the seismic imaging travel time data based on the pre-established upwind difference format data, the seismic velocity structure model, and the eikonal equation includes:

[0019] The eikonal equation is solved according to the basic data of the target work area, the upwind difference format data and the earthquake velocity structure model to determine the first arrival time of the earthquake.

[0020] In one embodiment, solving the eikonal equation based on the basic data of the target work area, the upwind difference format data, and the seismic velocity structure model to determine the earthquake first arrival time includes:

[0021] Determining the seismic ray emission direction of the seismic imaging data at the station according to the upwind difference format data;

[0022] performing a projection transformation operation on the seismic imaging data according to the direction of the seismic ray emission;

[0023] The long-short window ratio method is used to determine the ratio method arrival time based on the seismic data after projection transformation;

[0024] The backtracking travel time search algorithm is used to solve the eikonal equation according to the basic data, the earthquake velocity structure model and the ratio method arrival time to determine the earthquake first arrival time.

[0025] In one embodiment, the method for processing seismic imaging travel time data further includes: screening the seismic imaging data according to magnitude and time characterization parameters of the seismic imaging data.

[0026] In a second aspect, the present invention provides a seismic imaging travel time data processing device, comprising:

[0027] A structural model building module, configured to build a seismic velocity structural model of a target work area based on seismic imaging data of the target work area;

[0028] An eikonal equation building module, used for building the eikonal equation of the seismic imaging data;

[0029] The first arrival time determination module is used to determine the earthquake first arrival time of the seismic imaging travel time data based on the seismic velocity structure model, the eikonal equation and pre-established upwind difference format data.

[0030] In one embodiment, the eikonal equation establishment module includes:

[0031] The eikonal equation establishing unit is used to establish the eikonal equation according to the spatial travel time field, spatial gradient operator, medium slowness of the seismic imaging data and basic data of the target work area.

[0032] In one embodiment, the seismic imaging travel time data processing device further includes:

[0033] A format data establishment module is used to establish the upwind differential format data. The format data establishment module includes:

[0034] An operator determination unit, configured to determine an upwind difference operator based on the spatial travel time field of the seismic imaging data and its grid characterization parameters;

[0035] A format data establishing unit is used to establish the upwind difference format data according to the upwind difference operator.

[0036] In one embodiment, the first arrival time determination module includes:

[0037] The first arrival time determination unit is used to solve the eikonal equation based on the basic data of the target work area, the upwind difference format data and the seismic velocity structure model to determine the first arrival time of the earthquake.

[0038] In one embodiment, the first arrival time determining unit includes:

[0039] an emission direction determining unit, configured to determine the emission direction of the seismic ray of the seismic imaging data at the station according to the upwind differential format data;

[0040] A projection transformation unit, configured to perform a projection transformation operation on the seismic imaging data according to the direction of the seismic ray emission;

[0041] A ratio method arrival time determination unit is used to determine the ratio method arrival time based on the seismic data after projection transformation using the long-short window ratio method;

[0042] The earthquake first arrival time determination unit is used to solve the eikonal equation based on the basic data, the earthquake velocity structure model and the ratio method arrival time using a backtracking travel time search algorithm to determine the earthquake first arrival time.

[0043] In one embodiment, the seismic imaging travel time data processing device further includes:

[0044] The seismic data screening module is used to screen the seismic imaging data according to the magnitude and time characterization parameters of the seismic imaging data.

[0045] In a third aspect, the present invention provides an electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor implements the steps of a method for processing seismic imaging travel time data when executing the program.

[0046] In a fourth aspect, the present invention provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the seismic imaging travel-time data processing method.

[0047] As can be seen from the above description, the seismic imaging travel time data processing method and apparatus provided in the embodiments of the present invention first establish a seismic velocity structure model of the target work area based on the seismic imaging data of the target work area; then establish the Eikonal equation for the seismic imaging data; and finally determine the earthquake first arrival time of the seismic imaging travel time data based on the pre-established upwind difference format data, the seismic velocity structure model, and the Eikonal equation. The present invention addresses the technical pain points of low efficiency and poor precision in manual processing of seismic imaging travel time data in the prior art. By performing data screening, SAC file creation, theoretical travel time tagging, and automatic search of actual travel time, the present invention achieves fully automatic, efficient, and high-precision processing of seismic travel time data, providing a reliable technical guarantee for seismic travel time data processing for actual large-scale imaging problems. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0049] Figure 1This is a first structural diagram of a seismic imaging travel time data processing system according to an embodiment of the present application;

[0050] Figure 2 This is a second structural diagram of a seismic imaging travel time data processing system according to an embodiment of the present application;

[0051] Figure 3 Schematic diagram of the process of the seismic imaging travel time data processing method in an embodiment of the present invention Figure 1 ;

[0052] Figure 4 200 is a flow chart of step 200 in an embodiment of the present invention;

[0053] Figure 5 Schematic diagram of the process of the seismic imaging travel time data processing method in an embodiment of the present invention Figure 2 ;

[0054] Figure 6 4 is a flow chart of step 400 in an embodiment of the present invention;

[0055] Figure 7 300 is a flowchart of an embodiment of the present invention;

[0056] Figure 8 301 is a flow chart of step 301 in an embodiment of the present invention;

[0057] Figure 9 Schematic diagram of the process of the seismic imaging travel time data processing method in an embodiment of the present invention Figure 3 ;

[0058] Figure 10 Schematic diagram of the process of the seismic imaging travel time data processing method in an embodiment of the present invention Figure 4 ;

[0059] Figure 11 Schematic diagram of a flow chart of a method for processing travel-time data of seismic imaging in a specific application example of the present invention;

[0060] Figure 12 A regional velocity structure model (with a depth of 10 km) is established by screening based on a global model in a specific application example of the present invention;

[0061] Figure 13 A regional velocity structure model (with a depth of 20 km) established by screening on the basis of a global model in a specific application example of the present invention;

[0062] Figure 14 A regional velocity structure model (with a depth of 30 km) established by screening on the basis of a global model in a specific application example of the present invention;

[0063] Figure 15 A regional velocity structure model (with a depth of 40 km) is established by screening based on a global model in a specific application example of the present invention;

[0064] Figure 16 The event waveform diagram is formed by intercepting and processing the original continuous waveform in the specific application example of the present invention;

[0065] Figure 17 A schematic diagram of theoretical arrival time (N component of earthquake displacement) obtained by solving the eikonal equation in a specific application example of the present invention;

[0066] Figure 18 A schematic diagram of theoretical arrival time (E component of earthquake displacement) obtained by solving the Eikonal equation in a specific application example of the present invention;

[0067] Figure 19 A schematic diagram of theoretical arrival time (Z component of earthquake displacement) obtained by solving the Eikonal equation in a specific application example of the present invention;

[0068] Figure 20 A schematic diagram of the seismic wave displacement transformation along the seismic exit square projection and the ratio of the long-short time window in a specific application example of the present invention;

[0069] Figure 21 Schematic diagram of earthquake waveforms recorded by 35 broadband seismographs in Tengchong and its surrounding areas in a specific application example of the present invention;

[0070] Figure 22 1 is a structural block diagram of a seismic imaging travel time data processing device in an embodiment of the present invention;

[0071] Figure 23 Schematic diagram of the structure of an electronic device in an embodiment of the present invention. DETAILED DESCRIPTION

[0072] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention without making any creative efforts shall fall within the scope of protection of the present invention.

[0073] In one embodiment, the present application also provides a seismic imaging travel time data processing system, see Figure 1 The system can be a server A1, which can be connected to multiple stations (referring to observation points using various seismic instruments for earthquake observation) B1 in communication. The server A1 can also be connected to multiple databases in communication, or as shown in FIG. Figure 2 As shown, these databases can also be directly set up in server A1. Sensor B1 is used to measure teleseismic data in real time. After receiving the teleseismic data, server A1 performs preprocessing and seismic imaging travel time data processing on the teleseismic data, and then displays the processing results (earthquake first arrival time) to the user on client C1.

[0074] It is understood that the station B1 may be a pressure sensor, a seismic data receiver, a flow sensor, a displacement sensor, etc., and the client C1 may include a smart phone, a tablet electronic device, a network set-top box, a portable computer, a desktop computer, a personal digital assistant (PDA), a vehicle-mounted device, a smart wearable device, etc. Among them, the smart wearable device may include smart glasses, smart watches, smart bracelets, etc.

[0075] In practical applications, the pre-processing and seismic imaging travel time data processing can be performed on the server A1 side as described above, that is, Figure 1 or Figure 2 The architecture shown can also be used in which all operations are completed in the client C1 device. The specific selection can be based on the processing capabilities of the client device and the limitations of the user's usage scenario. This application does not limit this. If all operations are completed in the client device, the client device may also include a processor for performing operations such as pre-processing of teleseismic data and processing of seismic imaging travel time data.

[0076] The client C1 device described above may include a communication module (i.e., a communication unit) that can establish a communication connection with a remote server to implement data transmission with the server. The server may include a server for preprocessing and seismic imaging travel time data processing. Other implementation scenarios may also include an intermediate platform server, such as a server on a third-party server platform that is communicatively linked to a prediction server for preprocessing and seismic imaging travel time data processing. The server may include a single computer device, a server cluster consisting of multiple servers, or a server structure of a distributed device.

[0077] The server and the client device may communicate using any suitable network protocol, including network protocols that have not yet been developed as of the filing date of this application. Examples of network protocols include TCP / IP, UDP / IP, HTTP, and HTTPS. Examples of network protocols include RPC (Remote Procedure Call Protocol) and REST (Representational State Transfer) that are used on top of the aforementioned protocols.

[0078] It should be noted that, in the absence of conflict, the embodiments and features of the embodiments in this application can be combined with each other. The present application will be described in detail below with reference to the accompanying drawings and in combination with the embodiments.

[0079] In the prior art, specific seismic phases are often calibrated manually in travel time tomography data processing. This method of calibrating seismic data processing is inefficient and consumes a lot of human resources. In order to address the drawbacks of low efficiency in manual calibration of earthquake arrival times, those skilled in the art have also tried to automatically calibrate arrival times through computers, that is, to quickly calibrate the arrival times of seismic waves using the seismic long-short time window ratio method. However, the long-short time window ratio method has obvious defects, and its defects are mainly manifested in two aspects. First, the long-short time window ratio method has poor accuracy in calibrating travel times. Although the long-short time window ratio method has high accuracy in calibrating arrival times for synthetic data, it has poor accuracy when used for actual data. The main reason is that actual seismic data is more complex and contains more irregular signals. The presence of these irregular signals interferes with the accurate calibration of earthquake travel times by the long-short time window ratio method. Based on the problems that seismic data travel time calibration requires a lot of human resources and traditional automatic calibration is not accurate enough, an embodiment of the present invention provides a specific implementation method of a seismic imaging travel time data processing method, see. Figure 3 , the method specifically includes the following contents:

[0080] Step 100: Establishing a seismic velocity structure model of the target work area based on seismic imaging data of the target work area.

[0081] Specifically, the regional earthquake catalog and the distribution of seismic stations in the region are obtained according to the selected seismic imaging research scope (target work area), and the seismic velocity structure model of the target work area is established based on the global model (such as the AK135 model) and the crustal correction model (such as CRUST1.0).

[0082] Step 200: Establish an eikonal equation for the seismic imaging data.

[0083] Specifically, the eikonal equation is established according to formula (1):

[0084]

[0085] Where T is the space travel time field, is the spatial gradient operator, S is the medium slowness; r, θ, are the three direction coordinates of spherical coordinates, r represents the distance between any point in space and the center of the earth, θ represents the co-latitude, Indicates longitude.

[0086] It can be understood that the eikonal equation established according to formula (1) is more suitable for obtaining travel time data of a seismic data model that is heterogeneous and has a complex structure.

[0087] Step 300: Determine the earthquake first arrival time of the seismic imaging travel time data according to the seismic velocity structure model, the eikonal equation and pre-established upwind difference format data.

[0088] It should be noted that the upwind difference format data in step 300 satisfies the conservation of entropy. When step 300 is implemented, specifically, the eikonal equation is solved based on the upwind difference format data, the seismic velocity structure model and the basic data of the target work area to determine the time of the first arrival of the earthquake.

[0089] As can be seen from the above description, the method for processing seismic imaging travel time data provided by an embodiment of the present invention first quickly screens out seismic event information and initial velocity structure model information for seismic imaging based on the global earthquake catalog, station information, and velocity structure model. Secondly, preliminary data processing is performed on the original continuous waveform data based on the event information to obtain the event waveform. Then, based on the event information and model information, the upwind difference format is used to solve the eikonal equation to determine the theoretical arrival time. Finally, a backtracking search algorithm is used to search for the precise earthquake first arrival time near the theoretical arrival time, and the seismic travel time data processing is completed. The present invention can replace manual data processing, realizes the full automation of seismic travel time data processing, and has the advantages of high efficiency and excellent accuracy.

[0090] In one embodiment, see Figure 4 , step 200 specifically includes:

[0091] Step 201: Establishing the eikonal equation based on the spatial travel time field, spatial gradient operator, medium slowness of the seismic imaging data and basic data of the target work area.

[0092] The specific implementation method is shown in formula (1), which will not be repeated here.

[0093] In one embodiment, see Figure 5 , the seismic imaging travel time data processing method further includes:

[0094] Step 400: Create the upwind difference format data, see Figure 6 , step 400 includes:

[0095] Step 401: determining an upwind difference operator according to the spatial travel time field of the seismic imaging data and its grid characterization parameters;

[0096] Step 402: Establish the upwind difference format data according to the upwind difference operator.

[0097] Specifically, step 401 and step 402 may be implemented according to formula (2):

[0098]

[0099] Where D represents the upwind difference operator, and the subscripts r, θ and Respectively represent the effects on r, θ and direction, the subscript of the travel time field T represents the grid point number; δr, δθ and is the spatial grid in r, θ and The grid discrete interval of the direction; r i Represents the distance between the i-th grid point and the center of the earth in the r direction.

[0100] In one embodiment, see Figure 7 , step 300 includes:

[0101] Step 301: Solve the Eikonal equation based on the basic data of the target work area, the upwind difference format data and the seismic velocity structure model to determine the first arrival time of the earthquake.

[0102] Specifically, based on the earthquake velocity structure model, the theoretical arrival time is determined by solving the Eikonal equation using basic data and upwind difference format data. Next, a backtracking search algorithm is used to search for the accurate earthquake first arrival time near the theoretical arrival time, completing the earthquake travel time data processing.

[0103] In one embodiment, see Figure 8 , step 301 includes:

[0104] Step 3011: determining the seismic ray emission direction of the seismic imaging data at the station according to the upwind differential format data;

[0105] The travel time field calculated according to formula (2) can be further used to obtain the direction of seismic ray emission at the station. The direction of seismic ray emission can be expressed as:

[0106]

[0107] where (e±1, f±1, g±1) are the eight grid points around the station; and are the three basis vectors of spherical coordinates; and are parallel to r, θ and Direction vector.

[0108] Step 3012: performing a projection transformation operation on the seismic imaging data according to the direction of the seismic ray emission;

[0109] According to formula (3), the three components of earthquake (in and are the three unit vectors along latitude, longitude and vertical at the station respectively; U N 、U Eand U Z The displacements on three unit vectors) are projected onto the ray exit direction The seismic waveform after projection transformation can be expressed as:

[0110]

[0111] where U l is the displacement in the seismic exit direction. It can be seen from formula (4) that the three components of the seismic wave are projected onto the seismic exit direction, and this transformation can effectively enhance the amplitude of the seismic wave.

[0112] Step 3013: Use the long-short window ratio method to determine the ratio method arrival time according to the seismic data after projection transformation;

[0113] If the filtering range of the seismic data is [f1, f2], the average period of the seismic wave signal is:

[0114]

[0115] The formula for determining the first arrival time by the traditional long-short time window ratio method is as follows:

[0116]

[0117] where 0 < a < b are constants, generally a is taken as 0.1 and b is taken as 3.0, τ i is any moment of the seismic record, max is the operation of taking the maximum value, and t3 is the arrival time calibrated by the traditional long-short time window ratio method.

[0118] Step 3014: Use the backtracking travel time search algorithm to solve the eikonal equation according to the basic data, the seismic velocity structure model, and the ratio method arrival time to determine the seismic first arrival time.

[0119] The travel time obtained from formula (6) is delayed compared to the true travel time. The fundamental reason is that the energy accumulated by the seismic signal for a period of time after the starting point can make the ratio of the long-short time window ratio increase and form a peak. In order to make the first arrival time calibrated by the long-short time window ratio closer to the true arrival time, in this embodiment, the backtracking travel time search algorithm is used, that is, backtracking forward according to the travel time t3 obtained from the long-short time window ratio. When the value of the long-short time window ratio is the average value of the long-short time window ratio, the corresponding travel time at this place is considered the accurate arrival time. The backtracking travel time search algorithm is expressed by the formula as follows:

[0120]

[0121] where abs is the operation of taking the absolute value, and t4 is the first arrival time obtained by backtracking search.

[0122] In one embodiment, refer to Figure 9, the seismic imaging travel time data processing method further includes:

[0123] Step 500: Filter the seismic imaging data according to the magnitude and time characterization parameters of the seismic imaging data.

[0124] According to regional seismic imaging principles, earthquakes are considered point sources. Therefore, the finite rupture surface of an earthquake is not considered, and seismic energy is assumed to be generated at a single point in space. Under conditions of relatively low magnitude (M 2.0-5.5), earthquakes can be considered point sources. One of the screening criteria for seismic imaging travel-time data is to select available earthquake events based on magnitudes of M 2.0-5.5. A large earthquake is typically preceded by a series of closely spaced foreshocks, and similarly, a series of closely spaced aftershocks follows. Given the possibility of a large number of earthquakes occurring in a given period, a second screening criterion for seismic imaging travel-time data is that the time difference between adjacent earthquakes must be greater than 60 seconds to prevent interference between them.

[0125] In one embodiment, see Figure 10 , the seismic imaging travel time data processing method further includes:

[0126] Step 600: Preprocess the seismic imaging travel time data.

[0127] Specifically, the filtering range is selected according to the needs of the actual research problem, and then the seismic imaging travel time data are subjected to detrending, deaveraging and despiking processing.

[0128] To further illustrate this solution, the present invention takes the Tengchong volcanic area as an example to provide a specific application example of the seismic imaging travel time data processing method. The specific application example includes the following contents, see Figure 11 .

[0129] S1: Screening of seismic imaging travel time data.

[0130] First, the research scope is determined. A catalog of earthquakes suitable for seismic imaging is selected based on the global earthquake catalog and the time period of station deployment. Based on this, a regional seismic velocity structure reference model is selected and established using global models (such as the AK135 model) and crustal correction models (such as CRUST1.0).

[0131] The regional earthquake catalog automatically screened by the global earthquake catalog is shown in Table 1. In Table 1, the first column is the earthquake event number; the second to eighth columns are the year, month, day, hour, minute, second, and millisecond of the earthquake occurrence time; the ninth to eleventh columns are the latitude, longitude, and depth of the earthquake, respectively.

[0132] The regional initial velocity model established by the global model is as follows Figures 12 to 15 As shown ( Figures 12 to 15The velocity structures at depths of 10km, 20km, 30km and 40km respectively. Table 1 shows that regional earthquakes mainly occur in the middle and upper crust at shallow depths. Figures 12 to 15 It can be seen that the regional velocity model established based on the global velocity model shows good lateral heterogeneity in the middle and upper crust, but shows strong lateral heterogeneity in the lower crust and the top of the upper mantle.

[0133] Table 1

[0134]

[0135] S2: Create SAC file.

[0136] It is understandable that making earthquake data into SAC format facilitates the storage of earthquake data. According to the time of earthquake occurrence provided by the earthquake catalog, the earthquake waveform data t1 before and t2 after the earthquake occurrence (t1 and t2 can be 100 seconds) are intercepted from the continuous waveform, and the intercepted earthquake waveform data are made into a SAC file according to the SAC file format. According to the earthquake longitude, latitude, and depth coordinates provided by the earthquake catalog, the earthquake coordinates are written into the header information of the SAC file. According to the longitude, latitude, and elevation coordinates of the station provided by the array catalog, the station coordinates are written into the header information of the SAC file. Finally, according to the time when the earthquake occurred, the time of earthquake occurrence is written into the header information of the SAC file.

[0137] Specifically, we select an earthquake event that occurred at 15:40:01 on January 2, 2012, with longitude, latitude, and depth of 98.4198°, 25.1320°, and 0 km, respectively, to introduce the actual earthquake data travel time processing. The original waveform data recorded by the seismic station HS.MIZ is as follows: Figure 16 As shown. Figure 16 It can be seen that before the data is processed, it is difficult to see useful event waveform signals ( Figure 16 (a) to Figure 16 (c)), so the first arrival time cannot be directly calibrated by the original waveform data. The SAC file production module is used to process the original waveform data file. According to the earthquake catalog, a section of earthquake signal is intercepted from the original waveform data to produce a SAC file, and the earthquake coordinates, earthquake time information and station information are written into the header information of the SAC file. Then, the earthquake signal is band-pass filtered, and finally the earthquake data is subjected to detrending, deaveraging and detrending processing, and the production of the SAC file is completed. The earthquake data signal after the SAC file is produced is as follows Figure 16 (d) to Figure 16 As shown in (f), the earthquake event waveform can be clearly seen.

[0138] S3: Mark the theoretical travel time.

[0139] Based on the earthquake coordinates and station coordinates, as well as the established regional seismic velocity reference model, the Eikonal equation is solved to obtain the theoretical reference travel time t3, and the obtained reference travel time is written into the header information of the SAC file. To be applicable to the calculation of synthetic travel time in any non-uniform structural model, the present invention adopts an upwind difference format that satisfies the conservation of entropy to solve the Eikonal equation with variable coefficients. This numerical format can stably obtain travel time in any non-uniform model with complex structures, ensuring the accuracy of the theoretical travel time mark. For details, see formula (1).

[0140] S4: Automatically search for actual time.

[0141] The specific implementation process is shown in formulas (3) to (7), which will not be repeated here.

[0142] S5: Result verification.

[0143] The theoretical travel time of the first arrival time is obtained by numerically solving the eikonal equation based on the selected earthquake catalog and the established velocity structure model. For the earthquake and station discussed in this invention, the calculated theoretical arrival time is 14.22 seconds, and the actual manually identified travel time is 14.49 seconds, a difference of only 0.27 seconds. This shows that the present invention can provide a more accurate theoretical arrival time by solving the eikonal equation, which is convenient for further marking more accurate first arrival times on this basis. Figures 17 to 19 To solve the eikonal equation to obtain the theoretical arrival time, the solid line in the figure represents the theoretical arrival time and the dotted line represents the actual arrival time. Figures 17 to 19 They correspond to the N, E and Z components of earthquake displacement respectively.

[0144] The earthquake displacement is projected according to the travel time search module, and the transformed radial displacement is as follows: Figure 20 As shown in (a), ( Figure 20 (a) is the waveform of the three components of earthquake displacement after being projected along the earthquake's outgoing direction. Figure 20 (b) is Figure 20 The waveform in (a) is compared with the ratio distribution diagram obtained by the long-short time window ratio method. Figure 20 (a) and Figures 17 to 19 It can be found that Figure 20 Figure (a) shows the projection transformation of the seismic displacement along the earthquake's exit square and the value of the long-short time window ratio. The waveforms of the three components of the seismic displacement after projection along the earthquake's exit direction are shown. The four black lines in the figure represent the theoretical travel time, the actual travel time, the travel time obtained by the backtracking search algorithm, and the travel time obtained by the long-short time window ratio.

[0145] After the radial projection transformation, the amplitude of the seismic wave is significantly enhanced, which is convenient for the subsequent precise marking of travel time. Figure 20 The waveform in (a) is the ratio distribution obtained by the long-short time window ratio method. Figure 20 As shown in (b). Figure 20 As can be seen from (b), near the initial arrival time, the long-short time window ratio has a local maximum of 0.378, and the corresponding marker arrival time is 15.06 seconds. Using the backtracking search algorithm, that is, the long-short time window ratio search considering the average noise level, the marker travel time obtained is 14.61 seconds. The theoretical arrival time, actual arrival time, long-short time window ratio marker arrival time, and backtracking search algorithm marker arrival time are plotted simultaneously on Figure 20 In (a), it can be seen that the arrival time error of the long-short time window ratio marker is large, reaching 0.57 seconds. The fundamental reason for the large error is the delay of the long-short time window ratio. Figure 20 As shown in Figure (a), the backtracking search algorithm obtains an initial arrival time of 14.61 seconds, which is the closest to the actual arrival time, with an error of only 0.12 seconds. This comparison shows that the backtracking search algorithm overcomes the delay caused by the long-short window ratio by searching forward and backward while taking noise into account, achieving high travel time marking accuracy.

[0146] In order to further verify the accuracy of the backtracking search algorithm in calibrating the actual data travel time, 35 additional seismic stations around Tengchong were selected for algorithm testing. The displacement of the earthquake signal emission direction recorded by these 35 seismic stations is as follows: Figure 21 As shown. Figure 21 It can be seen that the signal-to-noise ratio of the seismic signal is widely distributed, with seismic data of high signal-to-noise ratio (such as stations X1.53087 and X1.53088) and seismic data of low signal-to-noise ratio (such as stations HS.RST and X1.53040). Selecting seismic data with different seismic signal-to-noise ratios can better illustrate the processing capability of the backtracking search algorithm of the present invention for seismic signals of different qualities. The test results are shown in Table 2. For comparison purposes, Table 2 shows the accuracy of the travel time marking of the backtracking search algorithm in addition to the accuracy of the long-short time window ratio method. It can be seen from Table 2 that the travel time marking accuracy of the backtracking search algorithm is much higher than that of the traditional long-short time window ratio method. The average travel time error of the 35 stations marked by the backtracking search algorithm is 0.17 seconds, and the average value of the long-short time window ratio method is 0.46 seconds; the standard deviation of the backtracking search algorithm is 0.18, and the standard deviation of the long-short time window ratio method is 0.19. The comparison of the mean and standard deviation of the errors further confirms that the backtracking search algorithm has a high travel time marking capability and can meet the needs of actual seismic data imaging.

[0147] Table 2

[0148]

[0149]

[0150] As can be seen from the above description, the method for processing seismic imaging travel time data provided by the specific application example of the present invention first obtains a regional earthquake catalog, the distribution of seismic stations within the region, and the regional initial underground velocity structure model based on the selected seismic imaging research scope. Secondly, based on the earthquake catalog, seismic waveform data at t1 before and t2 after the earthquake (t1 and t2 can be 100 seconds) are intercepted from the continuous waveform data and generated into a SAC format file. The longitude, latitude, and depth coordinates of the earthquake and station are then written into the header information of the SAC file, and the earthquake onset time is also written into the SAC header information. Thirdly, based on the earthquake catalog and station coordinates and the underground velocity structure model, the eikonal equation is numerically solved to obtain the theoretical arrival time t3 of the earthquake at the station, and the theoretical arrival time is written into the header information of the SAC file. Finally, a backtracking search algorithm is used to automatically search for the earthquake arrival time near t3, and the travel time obtained from the search is written into the SAC header information, completing the processing of the seismic imaging travel time data.

[0151] Based on the same inventive concept, the embodiments of the present application also provide a seismic imaging travel time data processing device, which can be used to implement the method described in the above embodiments, such as the following embodiments. Since the principle of solving the problem by the seismic imaging travel time data processing device is similar to that of the seismic imaging travel time data processing method, the implementation of the seismic imaging travel time data processing device can refer to the implementation of the seismic imaging travel time data processing method, and the repeated parts will not be repeated. As used below, the term "unit" or "module" can be a combination of software and / or hardware that implements a predetermined function. Although the system described in the following embodiments is preferably implemented in software, implementation by hardware, or a combination of software and hardware, is also possible and conceived.

[0152] The embodiment of the present invention provides a specific implementation of a seismic imaging travel time data processing device capable of implementing a seismic imaging travel time data processing method, see Figure 22 The seismic imaging travel time data processing device specifically includes the following contents:

[0153] A structural model building module 10 is used to build a seismic velocity structural model of the target work area based on the seismic imaging data of the target work area;

[0154] An eikonal equation building module 20, configured to build an eikonal equation for the seismic imaging data;

[0155] The first arrival time determination module 30 is used to determine the earthquake first arrival time of the seismic imaging travel time data according to the seismic velocity structure model, the eikonal equation and pre-established upwind difference format data.

[0156] In one embodiment, the eikonal equation establishment module includes:

[0157] The eikonal equation establishing unit is used to establish the eikonal equation according to the spatial travel time field, spatial gradient operator, medium slowness of the seismic imaging data and basic data of the target work area.

[0158] In one embodiment, the seismic imaging travel time data processing device further includes:

[0159] A format data establishment module is used to establish the upwind differential format data. The format data establishment module includes:

[0160] An operator determination unit, configured to determine an upwind difference operator based on the spatial travel time field of the seismic imaging data and its grid characterization parameters;

[0161] A format data establishing unit is used to establish the upwind difference format data according to the upwind difference operator.

[0162] In one embodiment, the first arrival time determination module includes:

[0163] The first arrival time determination unit is used to solve the eikonal equation based on the basic data of the target work area, the upwind difference format data and the seismic velocity structure model to determine the first arrival time of the earthquake.

[0164] In one embodiment, the root first arrival time determining unit includes:

[0165] an emission direction determining unit, configured to determine the emission direction of the seismic ray of the seismic imaging data at the station according to the upwind differential format data;

[0166] A projection transformation unit, configured to perform a projection transformation operation on the seismic imaging data according to the direction of the seismic ray emission;

[0167] A ratio method arrival time determination unit is used to determine the ratio method arrival time based on the seismic data after projection transformation using the long-short window ratio method;

[0168] The earthquake first arrival time determination unit is used to solve the eikonal equation based on the basic data, the earthquake velocity structure model and the ratio method arrival time using a backtracking travel time search algorithm to determine the earthquake first arrival time.

[0169] In one embodiment, the seismic imaging travel time data processing device further includes:

[0170] The seismic data screening module is used to screen the seismic imaging data according to the magnitude and time characterization parameters of the seismic imaging data.

[0171] As can be seen from the above description, the seismic imaging travel time data processing device provided by the embodiment of the present invention first establishes a seismic velocity structure model of the target work area based on the seismic imaging data of the target work area; then establishes the Eikonal equation for the seismic imaging data; and finally determines the earthquake first arrival time of the seismic imaging travel time data based on the pre-established upwind difference format data, the seismic velocity structure model, and the Eikonal equation. The present invention addresses the technical pain points of low efficiency and poor precision in the manual processing of seismic imaging travel time data in the prior art. By screening data, creating SAC files, marking theoretical travel times, and automatically searching actual travel times, the present invention achieves fully automatic, efficient, and high-precision processing of seismic travel time data, providing a reliable technical guarantee for the processing of seismic travel time data for actual large-scale imaging problems.

[0172] The embodiments of the present application also provide a specific implementation of an electronic device capable of implementing all steps in the seismic imaging travel time data processing method in the above embodiment, see Figure 23 , electronic equipment specifically includes the following:

[0173] Processor 1201, memory 1202, communications interface 1203, and bus 1204;

[0174] The processor 1201 , the memory 1202 , and the communication interface 1203 communicate with each other via the bus 1204 ; the communication interface 1203 is used to implement information transmission between server-side devices, sensors, client devices, and other related devices.

[0175] The processor 1201 is used to call the computer program in the memory 1202. When the processor executes the computer program, all steps of the seismic imaging travel time data processing method in the above embodiment are implemented. For example, when the processor executes the computer program, the following steps are implemented:

[0176] Step 100: screening the seismic data according to pre-received time characterization data to determine waveform data of a target teleseismic event;

[0177] Step 200: Calculating relative travel times of waveform data between multiple stations based on a cross-correlation function of waveform data between the multiple stations;

[0178] Step 300: Calculate the relative travel time error based on the waveform data, the cross-correlation function, and the relative travel time.

[0179] The embodiments of the present application also provide a computer-readable storage medium capable of implementing all steps of the method for processing seismic imaging travel time data in the above embodiments. The computer-readable storage medium stores a computer program. When the computer program is executed by a processor, all steps of the method for processing seismic imaging travel time data in the above embodiments are implemented. For example, when the processor executes the computer program, the following steps are implemented:

[0180] Step 100: screening the seismic data according to pre-received time characterization data to determine waveform data of a target teleseismic event;

[0181] Step 200: Calculating relative travel times of waveform data between multiple stations based on a cross-correlation function of waveform data between the multiple stations;

[0182] Step 300: Calculate the relative travel time error based on the waveform data, the cross-correlation function, and the relative travel time.

[0183] The various embodiments in this specification are described in a progressive manner. Similar portions between the various embodiments can be referenced to each other. Each embodiment focuses on the differences between the other embodiments. In particular, the hardware + program embodiments are generally similar to the method embodiments, so their description is relatively simple. For relevant portions, refer to the description of the method embodiments.

[0184] The foregoing description of this specification describes specific embodiments. Other embodiments are within the scope of the appended claims. In some cases, the actions or steps recited in the claims can be performed in an order different from that described in the embodiments and still achieve the desired results. Furthermore, the processes depicted in the accompanying drawings do not necessarily require the specific order shown or the sequential order to achieve the desired results. In certain embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0185] It will be understood by those skilled in the art that embodiments of the present invention may be provided as methods, systems, or computer program products. Thus, the present invention may take the form of an entirely hardware embodiment, an entirely software embodiment, or an embodiment combining software and hardware. Furthermore, the present invention may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to magnetic disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0186] The present invention is described with reference to flowcharts and / or block diagrams of methods, devices (systems), and computer program products according to embodiments of the present invention. It should be understood that each process and / or block in the flowcharts and / or block diagrams, as well as combinations of processes and / or blocks in the flowcharts and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowcharts and / or block diagrams. Figure 1 a process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0187] These computer program instructions may also be stored in a computer readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 a process or multiple processes and / or boxes Figure 1 The function specified in one or more boxes.

[0188] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operational steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing the instructions executed on the computer or other programmable device for implementing the process. Figure 1 a process or multiple processes and / or boxes Figure 1 A step that specifies a function in one or more boxes.

[0189] Specific embodiments are used in the present invention to illustrate the principles and implementation methods of the present invention. The description of the above embodiments is only used to help understand the method of the present invention and its core ideas. At the same time, for those skilled in the art, according to the ideas of the present invention, there may be changes in the specific implementation methods and application scopes. In summary, the contents of this specification should not be understood as limiting the present invention.

Claims

1. A method for processing seismic imaging travel time data, characterized in that: include: Establishing a seismic velocity structure model of the target work area based on seismic imaging data of the target work area; Establishing an eikonal equation for the seismic imaging data; Determining the earthquake first arrival time of the seismic imaging travel time data according to the seismic velocity structure model, the eikonal equation and pre-established upwind difference format data; The eikonal equation for establishing the seismic imaging data comprises: The eikonal equation is established based on the spatial travel time field, spatial gradient operator, medium slowness of the seismic imaging data and basic data of the target work area, specifically: According to formula (1), the eikonal equation is established: Where T is the space travel time field, is the spatial gradient operator, S is the medium slowness; r, θ, are the three direction coordinates of spherical coordinates, r represents the distance between any point in space and the center of the earth, θ represents the co-latitude, Indicates longitude; The method for establishing the upwind difference format data includes: Determining an upwind difference operator according to the spatial travel time field of the seismic imaging data and its grid characterization parameters; The upwind difference format data is established according to the upwind difference operator, specifically according to formula (2): Where D represents the upwind difference operator, and the subscripts r, θ and Respectively represent the effects on r, θ and direction, the subscript of the travel time field T represents the grid point number; δr, δθ and is the spatial grid in r, θ and The grid discrete interval of the direction; r i Indicates the distance between the i-th grid point and the center of the earth in the r direction; Determining the earthquake first arrival time of the seismic imaging travel time data according to the seismic velocity structure model, the eikonal equation, and pre-established upwind difference format data includes: Solving the eikonal equation according to the basic data of the target work area, the upwind difference format data and the seismic velocity structure model to determine the first arrival time of the earthquake; Solving the eikonal equation according to the basic data of the target work area, the upwind difference format data, and the seismic velocity structure model to determine the earthquake first arrival time includes: The seismic ray emission direction of the seismic imaging data at the station is determined according to the upwind differential format data, specifically: The travel time field calculated according to formula (2) can be further used to obtain the direction of seismic ray emission at the station. The direction of seismic ray emission can be expressed as: where (e±1, f±1, g±1) are the eight grid points around the station; and are the three basis vectors of spherical coordinates; and are parallel to r, θ and direction vector; According to the direction of the seismic ray emission, a projection transformation operation is performed on the seismic imaging data, specifically: The three components of earthquake are converted into in and are the three unit vectors along latitude, longitude and vertical at the station respectively; U N 、U E and U Z is the displacement of three unit vectors; projected to the direction of ray emission The seismic waveform after projection transformation is expressed as: Among them, U l is the displacement in the direction of earthquake emission; The long-short window ratio method is used to determine the ratio method arrival time based on the seismic data after projection transformation. Specifically: If the filtering range of seismic data is [f1, f2], the average period of the seismic wave signal is: The formula for determining the first arrival time using the traditional long-short time window ratio method is as follows: where \(0 < a < b\) are constants, generally \(a = 0.1\) and \(b = 3.0\), \(\tau\) i is any moment of the seismic record, max is the operation of taking the maximum value, and \(t_3\) is the arrival time calibrated by the traditional short-long time window ratio method; The eikonal equation is solved using a backtracking traveltime search algorithm based on the basic data, the earthquake velocity structure model, and the ratio method arrival time to determine the earthquake first arrival time. Specifically: The travel time t3 obtained based on the long-short time window ratio is traced back. When the value of the long-short time window ratio is the average long-short time window ratio, the corresponding travel time is considered to be accurate. The backtracking travel time search algorithm is expressed as follows: Where abs is the absolute value operation, and t4 is the first arrival time obtained by backtracking search.

2. The method for processing seismic imaging travel time data according to claim 1, wherein: Also includes: The seismic imaging data is filtered according to the magnitude and time characterization parameters of the seismic imaging data.

3. A seismic imaging travel time data processing device, characterized in that: include: A structural model building module, configured to build a seismic velocity structural model of a target work area based on seismic imaging data of the target work area; An eikonal equation building module, used for building the eikonal equation of the seismic imaging data; a first arrival time determination module, configured to determine the earthquake first arrival time of the seismic imaging travel time data based on the seismic velocity structure model, the eikonal equation, and pre-established upwind difference format data; The eikonal equation building module includes: The eikonal equation establishing unit is used to establish the eikonal equation based on the spatial travel time field, spatial gradient operator, medium slowness of the seismic imaging data and basic data of the target work area, specifically: According to formula (1), the eikonal equation is established: Where T is the space travel time field, is the spatial gradient operator, S is the medium slowness; r, θ, are the three direction coordinates of spherical coordinates, r represents the distance between any point in space and the center of the earth, θ represents the co-latitude, Indicates longitude; The method for establishing the upwind difference format data includes: Determining an upwind difference operator according to the spatial travel time field of the seismic imaging data and its grid characterization parameters; The upwind difference format data is established according to the upwind difference operator, specifically according to formula (2): Where D represents the upwind difference operator, and the subscripts r, θ and Respectively represent the effects on r, θ and direction, the subscript of the travel time field T represents the grid point number; δr, δθ and is the spatial grid in r, θ and The grid discrete interval of the direction; r i Indicates the distance between the i-th grid point and the center of the earth in the r direction; Determining the earthquake first arrival time of the seismic imaging travel time data according to the seismic velocity structure model, the eikonal equation, and pre-established upwind difference format data includes: Solving the eikonal equation according to the basic data of the target work area, the upwind difference format data and the seismic velocity structure model to determine the first arrival time of the earthquake; Solving the eikonal equation according to the basic data of the target work area, the upwind difference format data, and the seismic velocity structure model to determine the earthquake first arrival time includes: The seismic ray emission direction of the seismic imaging data at the station is determined according to the upwind differential format data, specifically: The travel time field calculated according to formula (2) can be further used to obtain the direction of seismic ray emission at the station. The direction of seismic ray emission can be expressed as: where (e±1, f±1, g±1) are the eight grid points around the station; and are the three basis vectors of spherical coordinates; and are parallel to r, θ and direction vector; According to the direction of the seismic ray emission, a projection transformation operation is performed on the seismic imaging data, specifically: The three components of earthquake are converted into in and are the three unit vectors along latitude, longitude and vertical at the station respectively; U N 、U E and U Z is the displacement of three unit vectors; projected to the direction of ray emission The seismic waveform after projection transformation is expressed as: Among them, U l is the displacement in the direction of earthquake emission; The long-short window ratio method is used to determine the ratio method arrival time based on the seismic data after projection transformation. Specifically: If the filtering range of seismic data is [f1, f2], the average period of the seismic wave signal is: The formula for determining the first arrival time using the traditional long-short time window ratio method is as follows: where \(0 < a < b\) are constants, generally \(a = 0.1\) and \(b = 3.0\), \(\tau\) i is any moment of the seismic record, max is the operation of taking the maximum value, and \(t_3\) is the arrival time calibrated by the traditional long - short time window ratio method; The eikonal equation is solved using a backtracking traveltime search algorithm based on the basic data, the earthquake velocity structure model, and the ratio method arrival time to determine the earthquake first arrival time. Specifically: The travel time t3 obtained based on the long-short time window ratio is traced back. When the value of the long-short time window ratio is the average long-short time window ratio, the corresponding travel time is considered to be accurate. The backtracking travel time search algorithm is expressed as follows: Where abs is the absolute value operation, and t4 is the first arrival time obtained by backtracking search.

4. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein: When the processor executes the program, the steps of the seismic imaging travel time data processing method according to any one of claims 1 to 2 are implemented.

5. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps of the seismic imaging travel time data processing method according to any one of claims 1 to 2 are implemented.

Images