CN114398780A

CN114398780A - Seismic imaging travel data processing method and device

Info

Publication number: CN114398780A
Application number: CN202210024491.8A
Authority: CN
Inventors: 刘少林; 杨顶辉; 汪文帅; 申文豪
Original assignee: Tsinghua University; National Institute of Natural Hazards
Current assignee: Tsinghua University; National Institute of Natural Hazards
Priority date: 2022-01-10
Filing date: 2022-01-10
Publication date: 2022-04-26

Abstract

The invention provides a method and a device for processing seismic imaging travel data, wherein the method for processing the seismic imaging travel data comprises the following steps: establishing a seismic velocity structure model of a target work area according to seismic imaging data of the target work area; establishing an equation of a path function of the seismic imaging data; and determining the seismic first arrival time of the seismic imaging travel time data according to the seismic velocity structure model, the equation of the equation and pre-established windward difference format data. The method and the device for processing the seismic imaging travel data can replace manual processing of the seismic imaging travel data, realize full automation of seismic imaging travel data processing, and have the advantages of high efficiency and high precision.

Description

Seismic imaging travel data processing method and device

Technical Field

The invention relates to the technical field of seismic data processing, in particular to a method and a device for processing seismic imaging travel time data.

Background

In the prior art, the exploration of the internal structure of the earth using a geophysical field is considered to be the most effective means for studying the internal structure of the earth. In particular, the internal structure of the earth can be described by seismic imaging techniques, it being understood that seismic imaging is the imaging of the internal structure of the earth with elastic wavefields excited by seismic. Seismic imaging can image parameters such as velocity structures, anisotropic structures, attenuation structures, and interface structures of the subsurface medium. Today, popular seismic imaging methods can be largely divided into two broad categories. One is wave equation based adjoint imaging; another type is ray-based travel time imaging.

For the first kind of methods, the adjoint imaging method based on the wave equation performs waveform cross correlation through a seismic source forward-propagating wave field and a station reverse-time backward-propagating wave field, and further constructs a Free t derivative of an error functional about model parameters, and after the Free t derivative is obtained, a numerical optimization method can be used, such as: a Steepest Descent Method (Steepest Description Method), a Conjugate Gradient Method (Conjugate Gradient Method), a Newton's Method (Newton's Method), a BFGS Method and the like, and the target model is gradually obtained by iterative optimization. Adjoint imaging based on the wave equation is used for constructing the Frechet derivative by solving a wave method, and the method has higher imaging resolution (the maximum resolution can reach half of the shortest wavelength) because the method contains the kinematic (such as travel time) and dynamic information (such as amplitude and finite frequency) information of seismic waves. However, numerical wave equations are required for either numerical solution of the source wavefield or the station back-propagation wavefield, and more computational resources are required for numerical solution of the seismic wave equations. Therefore, the cost of imaging is high with the wave equation, but the actual imaging efficiency is not high.

For the second method, the ray time-lapse imaging has higher imaging efficiency. At high frequency approximations, seismic wave propagation may be considered to be along ray propagation. Ray travel time imaging is mainly carried out through the sensitivity of ray travel time to a model to construct a travel time functional derivative relative to a model parameter Frechet. Calculating the seismic ray path is a key link of ray travel time imaging. There are two main methods for calculating ray paths in the industry, one is a pseudo-bending method, and the other is an equation solving method of an equation. The pseudo-bending method is characterized in that disturbance is given to a ray path between a seismic source and a station through connection according to a certain mode, the ray path is continuously corrected, and a finally calculated ray path is obtained. The pseudo-bending method involves solving a spatial first derivative of the model velocity, so that the model is required to be smooth, but the actual earth medium often has discontinuous discontinuities, and the pseudo-bending method is difficult to be applied to the seismic imaging problem in the actual complex medium. The equation solving method based on the equation of the function considers the influence of the heterogeneity of the medium on travel time, can be suitable for any heterogeneous medium model, and has high model adaptability. The numerical solution function equation has small calculation amount and does not need to be solved by a large-scale parallel computer. The imaging of the equation can reach better balance between the imaging precision and the cost.

When the equation of the equation is used for imaging, the inverted data object is a travel time residual, and the travel time residual can be obtained by the difference between the observed travel time and the synthesized travel time. The synthetic travel time is convenient to obtain, and after the travel time on the space grid point is obtained by solving the equation of the equation function, the synthetic travel time at the station can be obtained in a space interpolation mode. However, the acquisition of the observation travel time needs to be manually calibrated in actual data, and the precision of the travel time calibration in the actual data directly influences the imaging precision. Usually the calibration of the observed travel time is manually recognized by a trained person. Although the manual identification is high in precision, the imaging result is reliable. However, manual identification requires a large consumption of human resources, and the speed of manual identification is slow. The number of earthquake travel times needing to be calibrated in current earthquake imaging reaches the level of one hundred thousand or even millions, and the requirement of current earthquake mass data processing is not met only by manual identification.

In summary, there is a need for an automatic processing method for seismic travel data to improve the data processing capability in the field of seismic imaging.

Disclosure of Invention

Aiming at the problems in the prior art, the method and the device for processing the seismic imaging travel data can replace manual processing of the seismic imaging travel data, realize full automation of seismic imaging travel data processing, and have the advantages of high efficiency and high precision.

In order to solve the technical problems, the invention provides the following technical scheme:

in a first aspect, the present invention provides a seismic imaging travel time data processing method, including:

establishing a seismic velocity structure model of a target work area according to seismic imaging data of the target work area;

establishing an equation of a path function of the seismic imaging data;

and determining the seismic first arrival time of the seismic imaging travel time data according to the seismic velocity structure model, the equation of the equation and pre-established windward difference format data.

In one embodiment, the establishing the equation of the seismic imaging data comprises:

and establishing the equation of the program function according to the space travel time field, the space gradient operator, the medium slowness and the basic data of the target work area of the seismic imaging data.

In one embodiment, the method for establishing the windward difference format data comprises the following steps:

determining an upwind difference operator according to the space travel time field of the seismic imaging data and grid representation parameters of the space travel time field;

and establishing the windward difference format data according to the windward difference operator.

In an embodiment, the determining the seismic first arrival time of the seismic imaging travel time data according to the pre-established windward difference format data, the seismic velocity structure model and the equation of the equation function includes:

and solving the equation according to the basic data of the target work area, the windward difference format data and the seismic velocity structural model to determine the first arrival time of the earthquake.

In an embodiment, the solving the equation of the equation according to the basic data of the target work area, the windward difference format data and the seismic velocity structural model to determine the first arrival time of the earthquake includes:

determining the seismic ray emergence direction of the seismic imaging data at the station according to the windward difference format data;

performing projection transformation operation on the seismic imaging data according to the emergent direction of the seismic ray;

determining arrival time of a ratio method according to the seismic data after projection transformation by using a long-short window ratio method;

and solving the equation of the function according to the basic data, the seismic velocity structure model and the arrival time of the ratio method by utilizing a backtracking travel time search algorithm to determine the first arrival time of the earthquake.

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

In a second aspect, the present invention provides a seismic imaging travel time data processing apparatus, comprising:

the structural model establishing module is used for establishing a seismic velocity structural model of the target work area according to the seismic imaging data of the target work area;

the engineering function equation establishing module is used for establishing an engineering function equation of the seismic imaging data;

and the first arrival time determining module is used for determining the seismic first arrival time of the seismic imaging travel time data according to the seismic velocity structure model, the equation of the equation and pre-established windward difference format data.

In one embodiment, the equation establishing module comprises:

and the engineering function equation establishing unit is used for establishing the engineering function equation according to the space travel time field, the space gradient operator, the medium slowness and the basic data of the target work area of the seismic imaging data.

In one embodiment, the seismic imaging travel time data processing apparatus further includes:

the format data establishing module is used for establishing the windward difference format data, and comprises:

the operator determining unit is used for determining an upwind difference operator according to the space travel time field of the seismic imaging data and the grid representation parameters of the space travel time field;

and the format data establishing unit is used for establishing the windward difference format data according to the windward difference operator.

In one embodiment, the first arrival time determining module includes:

and the first arrival time determining unit is used for solving the equation of the equation according to the basic data of the target work area, the windward difference format data and the seismic velocity structural model so as to determine the seismic first arrival time.

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

the emergent direction determining unit is used for determining the emergent direction of the seismic ray of the seismic imaging data at the station according to the windward difference format data;

the projection transformation unit is used for performing projection transformation operation on the seismic imaging data according to the emergence direction of the seismic rays;

the specific value method arrival time determining unit is used for determining the specific value method arrival time according to the seismic data after the projection transformation by using a long-short window specific value method;

and the earthquake first arrival time determining unit is used for solving the equation of the equation according to the basic data, the earthquake speed structure model and the arrival time of the ratio method by utilizing a backtracking travel time searching algorithm so as to determine the earthquake first arrival time.

and the seismic data screening module is used for screening the seismic imaging data according to the magnitude and time characterization parameters of the seismic imaging data.

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

In a fourth aspect, the present invention provides a computer readable storage medium having stored thereon a computer program which, when executed by a processor, carries out the steps of a method of processing seismic imaging walk time data.

As can be seen from the above description, the method and apparatus for processing seismic imaging travel data according to the embodiments of the present invention first establish a seismic velocity structural model of a target work area according to seismic imaging data of the target work area; then establishing a path function equation of the seismic imaging data; and finally, determining the earthquake first arrival time of the earthquake imaging travel time data according to the pre-established windward difference format data, the earthquake speed structure model and the equation of the equation. Aiming at the technical pain points of low efficiency, poor precision and the like of manual processing of earthquake imaging travel data in the prior art, the invention realizes full-automatic high-efficiency and high-precision processing of the earthquake travel data by screening the data, making SAC files, marking theoretical travel time and automatically searching actual travel time, and provides reliable technical support for processing the earthquake travel data in the practical large-scale imaging problem.

Drawings

In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below, and it is obvious that the drawings in the following description are some embodiments of the present invention, and for those skilled in the art, other drawings can be obtained according to these drawings without creative efforts.

FIG. 1 is a first schematic diagram of a seismic imaging walktime data processing system according to an embodiment of the present disclosure;

FIG. 2 is a second schematic diagram of a seismic imaging walktime data processing system according to an embodiment of the present disclosure;

FIG. 3 is a first flowchart illustrating a seismic imaging travel time data processing method according to an embodiment of the present invention;

FIG. 4 is a flowchart of step 200 in an embodiment of the present invention;

FIG. 5 is a schematic flow chart of a seismic imaging travel time data processing method according to an embodiment of the invention;

FIG. 6 is a flowchart illustrating a step 400 according to an embodiment of the present invention;

FIG. 7 is a flowchart illustrating step 300 according to an embodiment of the present invention;

FIG. 8 is a flowchart illustrating step 301 according to an embodiment of the present invention;

FIG. 9 is a third schematic flow chart of a seismic imaging travel time data processing method according to an embodiment of the invention;

FIG. 10 is a fourth schematic flow chart of a seismic imaging travel time data processing method according to an embodiment of the present invention;

FIG. 11 is a flow chart of a seismic imaging travel time data processing method in an embodiment of the invention;

FIG. 12 is a regional velocity structure model (depth 10KM) built on the basis of a global model by screening in an embodiment of the present invention;

FIG. 13 is a regional velocity structure model (depth of 20KM) built on the basis of global model by screening in the embodiment of the present invention;

FIG. 14 is a regional velocity structure model (depth 30KM) built on the basis of a global model by screening in an embodiment of the present invention;

FIG. 15 is a regional velocity structure model (with a depth of 40KM) built on the basis of a global model by screening in an embodiment of the present invention;

FIG. 16 is a waveform diagram of an event generated by intercepting and processing an original continuous waveform in an embodiment of the present invention;

FIG. 17 is a schematic diagram of theoretical arrival time (seismic displacement N component) obtained by solving equation of function in a specific application example of the present invention;

FIG. 18 is a schematic diagram of theoretical arrival time (seismic displacement E component) obtained by solving an equation of an equation in a specific application example of the present invention;

FIG. 19 is a schematic diagram of theoretical arrival time (seismic displacement Z component) obtained by solving equation of function in an embodiment of the invention;

FIG. 20 is a diagram illustrating the transformation of the seismic wave displacement along the square projection of the seismic event and the ratio of the long and short windows in an embodiment of the present invention;

FIG. 21 is a schematic diagram of the seismic event waveforms recorded by soaks and 35 peripheral broadband seismographs in an embodiment of the invention;

FIG. 22 is a block diagram of a data processing apparatus for seismic imaging travel time in an embodiment of the invention;

fig. 23 is a schematic structural diagram of an electronic device in an embodiment of the invention.

Detailed Description

In order to make the objects, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are some, but not all, embodiments of the present invention. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.

In one embodiment, the present application further provides a seismic imaging travel data processing system, which may be a server a1, see fig. 1, where the server a1 may be communicatively connected to a plurality of stations (referring to observation points where seismic observation is performed by various seismic instruments) B1, the server a1 may be communicatively connected to a plurality of databases, respectively, or as shown in fig. 2, the databases may be disposed in the server a 1. Where sensor B1 is used to measure the teleseismic data in real time. After receiving the teleseismic data, the server a1 performs preprocessing, seismic imaging travel data processing, and the like on the teleseismic data, and displays a processing result (earthquake first arrival time) to the user at the client C1.

It is understood that station B1 may be a pressure sensor, seismic data receiver, flow sensor, displacement sensor, etc., and client C1 may include a smartphone, tablet, web set-top box, laptop, desktop, Personal Digital Assistant (PDA), in-vehicle device, smart wearable device, etc. Wherein, intelligence wearing equipment can include intelligent glasses, intelligent wrist-watch, intelligent bracelet etc..

In practical applications, the preprocessing and seismic imaging data processing parts can be performed on the side of the server a1 as described above, that is, the architecture shown in fig. 1 or fig. 2, and all operations can be completed in the client C1 device. The selection may be specifically performed according to the processing capability of the client device, the limitation of the user usage scenario, and the like. This is not a limitation of the present application. If all the operations are completed in the client device, the client device may further include a processor for performing operations such as preprocessing of the remote seismic data, data processing during seismic imaging, and the like.

The client C1 device may have a communication module (i.e., a communication unit) to communicate with a remote server for data transmission. The server may include a server on the preprocessing, seismic imaging travel data processing side, and in other implementations may include a server of an intermediate platform, such as a server of a third party server platform communicatively linked to a prediction server of the preprocessing, seismic imaging travel data processing. The server may comprise a single computer device, or may comprise a server cluster formed by a plurality of servers, or a server structure of a distributed device.

The server and client devices may communicate using any suitable network protocol, including network protocols not yet developed at the filing date of this application. The network protocols may include, for example, TCP/IP protocol, UDP/IP protocol, HTTP protocol, HTTPS protocol, and the like. Of course, the network Protocol may also include, for example, an RPC Protocol (Remote Procedure Call Protocol) used above the above Protocol, a REST Protocol (Representational State Transfer Protocol), and the like.

It should be noted that the embodiments and features of the embodiments in the present application may be combined with each other without conflict. The present application will be described in detail below with reference to the embodiments with reference to the attached drawings.

In the prior art, a specific seismic facies is usually calibrated in a manual mode in time-lapse tomography data processing, and the method for calibrating seismic data processing has low efficiency and consumes large human resources. Aiming at the defect of low efficiency of manually calibrating the arrival time of the earthquake, the technical personnel in the field also try to automatically calibrate the arrival time through a computer, namely, a short-time window comparison method of the earthquake is adopted to quickly calibrate the arrival time of the earthquake waves. However, the long-short time window ratio method has obvious defects, and the defects are mainly shown in two aspects. Firstly, the calibration travel time precision of the long-time window ratio method is poorer. Although the time calibration of the synthetic data long and short time window ratio method has higher precision, the time precision of the synthetic data long and short time window ratio method is poorer when the synthetic data long and short time window ratio method is used for actual data, the main reason is that the actual seismic data are more complex, the seismic data contain more irregular signals, and the accurate calibration of the long and short time window ratio method on seismic travel time is interfered by the existence of the irregular signals. Based on the problems that seismic data travel time calibration needs to consume a large amount of human resources and the traditional automatic calibration precision is insufficient, the embodiment of the invention provides a specific implementation mode of a seismic imaging travel time data processing method, and the specific implementation mode is shown in fig. 3, and the method specifically comprises the following steps:

step 100: and establishing a seismic velocity structural model of the target work area according to the seismic imaging data of the target work area.

Specifically, a regional seismic catalog and regional seismograph distribution are obtained according to a selected seismic imaging research range (target work area), and a seismic velocity structural model of the target work area is established according to a global model (such as an AK135 model) and a CRUST correction model (such as a CRUST 1.0).

Step 200: and establishing an equation of the seismic imaging data.

Specifically, an equation of the equation is established according to formula (1):

wherein T is a space travel time field,

is a spatial gradient operator, and S is a medium slowness; the sum of the r, the theta,

three directional coordinates, r, being spherical coordinatesShowing the distance between any point in space and the center of the earth, theta represents the weft remainders,

indicating longitude.

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

Step 300: and determining the seismic first arrival time of the seismic imaging travel time data according to the seismic velocity structure model, the equation of the equation and pre-established windward difference format data.

It should be noted that the windward difference format data in step 300 satisfies entropy conservation, and when step 300 is implemented, specifically: solving an equation of the equation according to the windward difference format data, the seismic velocity structure model and the basic data of the target work area to determine the first arrival time of the earthquake.

As can be seen from the above description, in the seismic imaging travel data processing method provided in the embodiment of the present invention, firstly, the seismic event information and the initial velocity structure model information for seismic imaging are quickly screened out according to the global seismic directory, the station information, and the velocity structure model. And secondly, performing primary data processing on the original continuous waveform data according to the event information to obtain an event waveform. And then, solving a function equation by adopting an upwind difference format according to the event information and the model information to determine a theoretical arrival time. And finally, searching the accurate earthquake first arrival time near the theoretical arrival time through a backtracking search algorithm, and finishing data processing during earthquake walking. The invention can replace manual data processing, realizes full automation of data processing during earthquake walking, and has the advantages of high efficiency and excellent precision.

In one embodiment, referring to fig. 4, step 200 specifically includes:

step 201: and establishing the equation of the program function according to the space travel time field, the space gradient operator, the medium slowness and the basic data of the target work area of the seismic imaging data.

The concrete implementation method is shown in formula (1), and the description is omitted.

In one embodiment, referring to fig. 5, the seismic imaging travel time data processing method further includes:

step 400: establishing the windward difference format data, see fig. 6, step 400 includes:

step 401: determining an upwind difference operator according to the space travel time field of the seismic imaging data and grid representation parameters of the space travel time field;

step 402: and establishing the windward difference format data according to the windward difference operator.

Specifically, step 401 and step 402 can be implemented according to formula (2):

where D denotes the windward difference operator, subscripts r, θ and

respectively represent the effects on r, theta and

direction, subscripts of the travel time field T represent grid point numbers; δ r, δ θ and

for the spatial grid at r, θ and

a grid of discrete intervals of direction; r is_iThe distance between the ith grid point in the r direction and the geocenter is shown.

In one embodiment, referring to fig. 7, step 300 includes:

step 301: and solving the equation according to the basic data of the target work area, the windward difference format data and the seismic velocity structural model to determine the first arrival time of the earthquake.

Specifically, according to the seismic velocity structural model, a theoretical arrival time is determined by solving an equation of a function of a equation by adopting basic data and windward difference format data. Then, accurate earthquake first arrival time is searched near the theoretical arrival time through a backtracking search algorithm, and data processing during earthquake walking is completed.

In one embodiment, referring to fig. 8, step 301 comprises:

step 3011: determining the seismic ray emitting direction of the seismic imaging data at the station according to the windward difference format data;

the travel time field calculated according to the formula (2) can further obtain the emergence direction of the seismic ray at the station, and the emergence direction of the seismic ray can be expressed as:

where (e ± 1, f ± 1, g ± 1) are 8 grid points around the station;

and

three basis vectors of the spherical coordinate are respectively;

and

are parallel to r, theta and

vector of direction.

Step 3012: performing projection transformation operation on the seismic imaging data according to the emergent direction of the seismic ray;

the earthquake three components can be obtained by the formula (3)

(wherein

And

three unit vectors along latitude, longitude and vertical at the station respectively; u shape_N、U_EAnd U_ZDisplacement in three unit vectors) to the ray exit direction

The projective transformed seismic waveform may be represented as:

wherein U is_lIs the displacement in the direction of the seismic emission. As can be seen from the formula (4), the three seismic components are projected to the seismic emergence direction, and the transformation can effectively strengthen the amplitude of the seismic wave.

Step 3013: determining arrival time of a ratio method according to the seismic data after projection transformation by using a long-short window ratio method;

if the filtering range of the seismic data is f₁,f₂]The average period of the seismic signals is:

the formula of determining the first arrival time by the traditional long-short time window ratio method is as follows:

wherein 0<a<b is a constant, generally a is 0.1, b is 3.0, tau_iAt any time of seismic recording, max is the maximum value operation, t₃The time is calibrated by a traditional long-short time window comparison method.

Step 3014: and solving the equation of the function according to the basic data, the seismic velocity structure model and the arrival time of the ratio method by utilizing a backtracking travel time search algorithm to determine the first arrival time of the earthquake.

The travel time obtained from equation (6) is delayed compared to the true travel time. The fundamental reason for this is that the energy accumulated by the seismic signal over a period of time after the takeoff point increases the ratio of the long-to-short window ratio and forms a peak. In order to make the long and short time windows closer to the real arrival time than the calibrated first arrival time, in this embodiment, a backtracking travel time search algorithm is used, that is, the travel time t obtained according to the long and short time window ratio₃And backtracking forwards, and when the value of the long-short time window ratio is the value of the average long-short time window ratio, considering the corresponding travel time at the position as the accurate arrival time. The backtracking walk time search algorithm is formulated as follows:

wherein abs is an absolute value operation, t₄The first arrival time is obtained for backtracking search.

In one embodiment, referring to fig. 9, the method for processing the seismic imaging travel time data further includes:

step 500: and screening the seismic imaging data according to the magnitude and time characterization parameters of the seismic imaging data.

According to the regional seismic imaging principle, the earthquake is considered as a point source, so the limited fracture surface of the earthquake is not considered, and only the seismic energy is considered to be excited at one point in space. Under the condition that the magnitude of the earthquake (M2.0-5.5) is small, the earthquake can be regarded as a point source. One of the screening conditions for seismic imaging travel time data is to screen for available seismic events according to magnitude M2.0-5.5. A series of close-spaced pre-earthquakes usually occur before a major earthquake, and also close-spaced aftershocks occur after a major earthquake. In order to avoid mutual interference between earthquakes, the second screening condition of the earthquake imaging travel time data is that the difference between the origin times of adjacent earthquakes is more than 60 seconds.

In one embodiment, referring to fig. 10, the seismic imaging travel time data processing method further includes:

step 600: and preprocessing the seismic imaging travel time data.

Specifically, the filtering range is selected according to the requirements of practical research problems, and then the seismic imaging travel time data is subjected to de-trending, de-averaging and de-pinch processing.

To further illustrate the present solution, the present invention provides a specific application example of the data processing method during seismic imaging walking, taking soar volcanic area as an example, and the specific application example includes the following contents, see fig. 11.

S1: and screening the seismic imaging travel time data.

Firstly, determining a research range, and screening the seismic catalog which can be used for seismic imaging according to the global seismic catalog and the station laying time period. And on the basis, screening and establishing a regional seismic velocity structure reference model according to a global model (such as an AK135 model) and an crustal rectification model (such as a CRUST 1.0).

The regional seismic catalog automatically screened by the global seismic catalog is shown in table 1, and in table 1, the first column is the seismic event number; the 2 nd to 8 th columns become the year, month, day, hour, minute, second and millisecond of the earthquake origin time; columns 9 to 11 are latitude, longitude and depth of the earthquake, respectively.

The area initial velocity models established by the global model are shown in fig. 12 to 15 (fig. 12 to 15 are velocity structures at depths of 10km, 20km, 30km, and 40km, respectively). It can be seen from table 1 that regional earthquakes mainly occur in the upper and middle crustae, with shallow depths. As can be seen from fig. 12 to 15, the regional velocity model established from the global velocity model exhibits better lateral non-uniformity in the middle and upper crust, but stronger lateral non-uniformity at the top of the lower crust and upper mantle.

TABLE 1

S2: a SAC file is made.

It will be appreciated that making the seismic data into SAC format facilitates storage of the seismic data. According to the origin time provided by the earthquake directory, the earthquake waveform data of the front t1 and the back t2 (100 seconds can be taken for t1 and t 2) of the origin time is intercepted from the continuous waveform, and the intercepted earthquake waveform data is made into an SAC file according to an SAC file format. And providing earthquake longitude, latitude and depth coordinates according to the earthquake directory, and writing the earthquake coordinates into header information of the SAC file. And writing the station coordinates into the header information of the SAC file according to the longitude, latitude and elevation coordinates of the station provided by the array directory. And finally, writing the earthquake starting time into the header information of the SAC file according to the earthquake occurrence time.

Specifically, seismic events that occur at 15 hours, 1 month, 2 days, 2012, 40 minutes and 1 second, and have a longitude, latitude and depth of 98.4198 °, 25.1320 ° and 0km, respectively, are selected to introduce actual seismic data travel time processing. Miz recorded raw waveform data for the seismic station hs.is shown in fig. 16. As can be seen from fig. 16, before data is not processed, it is difficult to see useful event waveform signals (fig. 16 (a) to 16 (c)), and thus the first arrival time cannot be directly scaled by the original waveform data. Processing an original waveform data file by adopting a SAC file manufacturing module, intercepting a section of seismic signals from the original waveform data according to an earthquake directory to manufacture a SAC file, writing seismic coordinates, earthquake-initiating time information and station information into header information of the SAC file, then performing band-pass filtering on the seismic signals, finally performing trend removing, averaging removing and trend removing processing on the seismic data, and finishing manufacturing the SAC file. As shown in fig. 16 (d) to 16 (f), the seismic data signals after the SAC file creation clearly show the waveform of the seismic event.

S3: the mark theory travel time.

Solving a function equation to obtain a theoretical reference travel time t according to the seismic coordinates and station coordinates and the established regional seismic velocity reference model₃And writing the obtained reference travel time into the header information of the SAC file. In order to be suitable for calculating the synthetic travel time in any non-uniform structure model, the invention adopts an upwind difference format which meets the entropy conservation to solve the equation of the variable coefficient, and the numerical format can stably obtain any non-uniform model and has travel time in a complex structure, thereby ensuring the theoretical travel timeThe accuracy of the travel time mark is shown in formula (1).

S4: the actual travel time is automatically searched.

The specific implementation process is shown in equations (3) to (7), which will not be described herein.

S5: and (5) verifying the result.

And obtaining the theoretical travel time of the first arrival time by solving an equation of the equation through numerical values according to the selected seismic catalogue and the established speed structure model. For the seismic and station discussed in this invention, the theoretical arrival time calculated is 14.22 seconds, and the actual manually identified travel time is 14.49 seconds, with a difference of only 0.27 seconds. Therefore, the invention can provide more accurate theoretical arrival time by solving the equation of the function of the equation, and is convenient for further marking more accurate initial arrival time on the basis. Fig. 17 to 19 show theoretical arrival times obtained by solving the equation of the equation function, in which the solid line represents the theoretical arrival time and the dotted line represents the actual arrival time. Fig. 17-19 correspond to the N, E and Z components of seismic displacement, respectively.

When the travel time search module projects the seismic displacement, the transformed radial displacement is as shown in fig. 20 (a), (fig. 20 (a) is a waveform obtained by projecting three components of the seismic displacement along the seismic emission direction, fig. 20 (b) is a ratio distribution diagram obtained by a long-short time window ratio method of the waveform in fig. 20 (a)) comparing fig. 20 (a) and fig. 17 to 19, and in fig. 20 (a), the seismic displacement is transformed along the seismic emission square projection and the values of the long-short time window ratio. The three components of the earthquake displacement are projected along the earthquake emergent direction to form a waveform, and four black straight lines in the diagram respectively represent the theoretical travel time, the real travel time, the travel time obtained by a backtracking search algorithm and the travel time obtained by a long-time window-to-short-time window ratio.

The amplitude of the seismic wave is obviously strengthened after the radial projection transformation, and the subsequent travel time accurate marking is facilitated. The waveform in fig. 20 (a) is shown in (b) of a ratio distribution chart 20 obtained by a long-short time window ratio method. As can be seen from fig. 20 (b), the long-short time window ratio has a local maximum value of 0.378 in the vicinity of the first arrival time, and the corresponding mark arrival time is 15.06 seconds. By using a backtracking search algorithm, namely searching by considering the ratio of the long time window to the short time window of the average noise level, the obtained mark travel time is 14.61 seconds. The theoretical arrival time, the actual arrival time, the long-short time window ratio marking arrival time and the backtracking search algorithm marking arrival time are simultaneously drawn in (a) of fig. 20, and it can be seen that the arrival time error of the long-short time window ratio marking is larger, the error reaches 0.57 second, and the root cause of the large error is the delay of the long-short time window ratio. As can also be seen from fig. 20 (a), the first arrival time obtained by the backtracking search algorithm is 14.61 seconds, which is the closest to the actual arrival time, and the error is only 0.12 seconds. The comparison result shows that the backtracking search algorithm carries out forward backtracking search under the condition of considering noise interference, overcomes the delay of the long and short window ratio and has higher travel time marking precision.

In order to further verify the accuracy of calibrating the actual data travel time by the backtracking search algorithm, 35 seismic stations around the flight-rush are additionally selected for algorithm inspection. The displacement of the seismic signal emitting directions recorded by the 35 seismic stations is shown in fig. 21. It can be seen from fig. 21 that the signal-to-noise ratio distribution of the seismic signals is broad, with seismic data with high signal-to-noise ratios (e.g., X1.53087 and X1.53088 stations) and seismic data with low signal-to-noise ratios (e.g., hs.rst and X1.53040 stations). The processing capability of the backtracking search algorithm of the invention on seismic signals with different qualities can be better illustrated by selecting seismic data with different seismic signal-to-noise ratios. The test results are shown in table 2. For comparison, table 2 shows the accuracy of the long-short time window comparison method in addition to the accuracy of the travel time marker of the backtracking search algorithm. As can be seen from table 2, the accuracy of the travel time marker of the backtracking search algorithm is much higher than that of the conventional long-short time window comparison method. The average value of the 35 station travel time errors marked by the backtracking search algorithm is 0.17 second, and the average value of the long-short time window comparison method is 0.46 second; 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 average value and the standard deviation of the error further proves that the backtracking search algorithm has higher travel time marking capability and can meet the requirement of actual seismic data imaging.

TABLE 2

As can be seen from the above description, the seismic imaging travel time data processing method provided in the specific application example of the present invention first obtains the regional seismic directory, the regional seismic station distribution, and the regional subsurface initial velocity structural model according to the selected seismic imaging research range. Secondly, seismic waveform data of a front t1 and a rear t2 (100 seconds for t1 and t 2) of the earthquake origin time are intercepted from the continuous waveform data according to an earthquake catalog and are made into SAC format files, and then longitude, latitude and depth coordinates of the earthquake and the station are written into header information of the SAC files, and the earthquake origin time is also written into SAC header information. Thirdly, according to the earthquake catalogue, the station coordinates and the underground velocity structure model, a mathematical solution function equation is carried out to obtain a theoretical arrival time t3 of the earthquake to the station, and the theoretical arrival time is written into the header information of the SAC file. And finally, automatically searching the earthquake arrival time near t3 according to a backtracking search algorithm, and writing the travel time obtained by searching into SAC header information so as to finish the data processing of the seismic imaging travel time.

Based on the same inventive concept, the embodiment of the present application further provides a seismic imaging travel time data processing apparatus, which can be used to implement the method described in the above embodiment, such as the following embodiments. Because the principle of solving the problems of the seismic imaging travel time data processing device is similar to 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 repeated parts are not described again. As used hereinafter, the term "unit" or "module" may be a combination of software and/or hardware that implements a predetermined function. While the system described in the embodiments below is preferably implemented in software, implementations in hardware, or a combination of software and hardware are also possible and contemplated.

The embodiment of the present invention provides a specific implementation of a seismic imaging travel data processing apparatus capable of implementing a seismic imaging travel data processing method, and referring to fig. 22, the seismic imaging travel data processing apparatus specifically includes the following contents:

the structural model establishing module 10 is used for establishing a seismic velocity structural model of the target work area according to the seismic imaging data of the target work area;

an equation establishing module 20, configured to establish an equation of the seismic imaging data;

and the first arrival time determining module 30 is configured to determine the seismic first arrival time of the seismic imaging travel time data according to the seismic velocity structural model, the equation of the equation and pre-established windward difference format data.

In one embodiment, the equation establishing module comprises:

In one embodiment, the first arrival time determining module includes:

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

the emergent direction determining unit is used for determining the emergent direction of the seismic rays of the seismic imaging data at the station according to the windward difference format data;

As can be seen from the above description, in the seismic imaging travel time data processing apparatus provided in the embodiment of the present invention, a seismic velocity structural model of a target work area is first established according to seismic imaging data of the target work area; then establishing a path function equation of the seismic imaging data; and finally, determining the earthquake first arrival time of the earthquake imaging travel time data according to the pre-established windward difference format data, the earthquake speed structure model and the equation of the equation. Aiming at the technical pain points of low efficiency, poor precision and the like of manual processing of earthquake imaging travel data in the prior art, the invention realizes full-automatic high-efficiency and high-precision processing of the earthquake travel data by screening the data, making SAC files, marking theoretical travel time and automatically searching actual travel time, and provides reliable technical support for processing the earthquake travel data in the practical large-scale imaging problem.

The embodiment of the present application further provides a specific implementation manner of an electronic device capable of implementing all steps in the seismic imaging travel time data processing method in the foregoing embodiment, and referring to fig. 23, the electronic device specifically includes the following contents:

a processor (processor)1201, a memory (memory)1202, a communication Interface 1203, and a bus 1204;

the processor 1201, the memory 1202 and the communication interface 1203 complete communication with each other through the bus 1204; the communication interface 1203 is used for implementing information transmission between related devices such as a server-side device, a sensor, a client device, and the like.

The processor 1201 is configured to call the computer program in the memory 1202, and the processor executes the computer program to implement all the steps in the seismic imaging walking data processing method in the above embodiments, for example, the processor executes the computer program to implement the following steps:

step 100: screening the seismic data according to the pre-received time characterization data of the seismic data to determine waveform data of a target remote earthquake event;

step 200: calculating relative travel time of waveform data among a plurality of stations according to a cross-correlation function of the waveform data among the stations;

step 300: and calculating the relative travel time error according to the waveform data, the cross-correlation function and the relative travel time.

Embodiments of the present application also provide a computer-readable storage medium capable of implementing all steps in the seismic imaging walking data processing method in the foregoing embodiments, where the computer-readable storage medium stores thereon a computer program, and when the computer program is executed by a processor, the computer program implements all steps of the seismic imaging walking data processing method in the foregoing embodiments, for example, when the processor executes the computer program, the processor implements the following steps:

The embodiments in the present specification are described in a progressive manner, and the same and similar parts among the embodiments are referred to each other, and each embodiment focuses on the differences from the other embodiments. In particular, for the hardware + program class embodiment, since it is substantially similar to the method embodiment, the description is simple, and the relevant points can be referred to the partial description of the method embodiment.

The foregoing description has been directed to specific embodiments of this disclosure. Other embodiments are within the scope of the following claims. In some cases, the actions or steps recited in the claims may be performed in a different order than in the embodiments and still achieve desirable results. In addition, the processes depicted in the accompanying figures do not necessarily require the particular order shown, or sequential order, to achieve desirable results. In some embodiments, multitasking and parallel processing may also be possible or may be advantageous.

As will be appreciated by one skilled in the art, embodiments of the present invention may be provided as a method, system, or computer program product. Accordingly, the present invention may take the form of an entirely hardware embodiment, an entirely software embodiment or an embodiment combining software and hardware aspects. Furthermore, the present invention may take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, and the like) having computer-usable program code embodied therein.

The present invention is described with reference to flowchart illustrations and/or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each flow and/or block of the flow diagrams and/or block diagrams, and combinations of flows and/or blocks in the flow diagrams and/or block diagrams, can be implemented by computer program instructions. These computer program instructions may be provided to a processor of a general purpose computer, special purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions specified in the flowchart flow or flows and/or block diagram block or blocks.

These computer program instructions may also be stored in a computer-readable memory that can direct a computer or other programmable data processing apparatus to function in a particular manner, such that the instructions stored in the computer-readable memory produce an article of manufacture including instruction means which implement the function specified in the flowchart flow or flows and/or block diagram block or blocks.

These computer program instructions may also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart flow or flows and/or block diagram block or blocks.

The principle and the implementation mode of the invention are explained by applying specific embodiments in the invention, and the description of the embodiments is only used for helping to understand the method and the core idea of the invention; meanwhile, for a person skilled in the art, according to the idea of the present invention, there may be variations in the specific embodiments and the application scope, and in summary, the content of the present specification should not be construed as a limitation to the present invention.

Claims

1. A seismic imaging travel time data processing method, comprising:

establishing an equation of a path function of the seismic imaging data;

2. The seismic imaging travel time data processing method of claim 1, wherein the establishing the equation of the path function of the seismic imaging data comprises:

3. The seismic imaging travel data processing method of claim 1, wherein the method of building the windward difference format data comprises:

4. The seismic imaging travel time data processing method of claim 1, wherein the determining the seismic first arrival time of the seismic imaging travel time data from the seismic velocity structure model, the equation of the function and pre-established windward difference format data comprises:

5. The seismic imaging travel time data processing method of claim 4, wherein solving the equation of the equation according to the basic data of the target work area, the windward difference format data and the seismic velocity structural model to determine the seismic first arrival time comprises:

6. The seismic imaging travel data processing method of claim 1, further comprising: and screening the seismic imaging data according to the magnitude and time characterization parameters of the seismic imaging data.

7. A seismic imaging travel time data processing apparatus, comprising:

8. The seismic imaging travel time data processing apparatus of claim 7 wherein the equation-of-function establishing module comprises:

9. An electronic device comprising a memory, a processor and a computer program stored on the memory and executable on the processor, wherein the processor when executing the program implements the steps of the seismic imaging walking data processing method of any of claims 1 to 6.

10. A computer-readable storage medium, on which a computer program is stored which, when being executed by a processor, carries out the steps of the seismic imaging walk time data processing method according to any one of claims 1 to 6.