Improved method, device, equipment and medium for inverting interval velocity by using first arrival waves

Abstract

The invention provides an improved method, device, equipment and medium for inverting interval velocity by using a first arrival wave, and the method comprises the steps: obtaining an actually measured wave field value, calculating the derivative of the actually measured wave field value, and picking up the first arrival time of a seismic record corresponding to the actually measured wave field value; determining an initial speed model; solving a wave field value calculated by using the current speed model, and calculating a derivative of the wave field value; calculating first arrival time from the excitation point and the receiving point to each grid point in the speed model; determining a propagation path of a first arrival wave; calculating to obtain a new speed model; whether the new speed model meets the precision is judged, and if the new speed model meets the precision, the current speed model is an inversion speed model; and if the precision is not met, repeating iteration to construct a new speed module until the new speed model meets the precision. According to the technical scheme, the precision of the inversion speed model can be improved to a large extent, the speed resolution is greatly improved, the anti-interference performance is high, and the interference of noise on seismic signals can be well suppressed.

Description

Technical Field

[0001] This application relates to the field of seismic data layer velocity inversion technology, and specifically to an improved method, apparatus, equipment and medium for inverting layer velocities using first arrival waves. Background Technology

[0002] In seismic exploration, accurately determining the velocity of subsurface strata is undoubtedly the most critical issue in the field of seismic data processing and interpretation. Based on the velocity requirements of different tasks, seismic exploration technicians generally classify velocities into several types, including average velocity, layer velocity, root-mean-square velocity, and stacking velocity. Among these, obtaining accurate layer velocities is key to constructing a subsurface velocity model. It can be said that if an accurate subsurface layer velocity model can be obtained, technicians can essentially determine the precise structural morphology of subsurface geological structures and the precise location of oil and gas reservoirs. Given the widespread application of surface seismic exploration technology in the field of seismic exploration, people usually use reflection wave information from 2D and 3D seismic data to determine the layer velocities of subsurface strata. The most commonly used method for determining layer velocities is to use horizontal stacking data and velocity spectrum data. Based on the principles of geometric seismology, assuming the subsurface strata are horizontally layered or horizontally inclined layered media, the calculation relationship between stacking velocity and layer velocity is used to recursively deduce the layer velocities of each subsurface layer from top to bottom. (Reference 1, Chen Chuanren & Zhou Xixiang) A precise method for inverting layer velocity. Geophysical and Geochemical Exploration Computation Technology, 1998, 3(21): 212-215) describes this method in detail. Since the actual underground geological conditions differ greatly from the assumptions given by this method in most cases, the layer velocity inverted by this method has a large error with the actual situation. Therefore, people usually use the layer velocity obtained by this method as the initial value, and then use other technical methods to improve the accuracy of the layer velocity based on the initial value. Many scientific and technological literatures have introduced this kind of technology. The technology introduced in Literature 2 (Wu Guochen, Wang Huazhong, Ma Zaitian. Constant velocity gradient ray tracing and two-dimensional layer velocity inversion. Petroleum Geophysical Exploration, 2003, 4(42): 434-440) is one of such technologies.

[0003] Generally speaking, current layer velocity inversion mainly utilizes reflected wave information from seismic data. This is because current seismic exploration primarily uses surface seismic exploration, and the seismic wave field information received at the surface can only be reflected wave information. However, with the gradual application of non-zero bias VSP seismic exploration technology and well-to-well tomography technology in the field of seismic exploration, how to utilize the transmitted wave field information to invert layer velocity has become an important issue in the field of seismic exploration technology. In the past period, researchers in the field of geophysical methods and technologies in exploration have done a lot of work in this area. Reference 3 (Y. Luo and T. Schuster, 1991, Wave-equation traveltime inversion, Geophysics, 56:645~653) proposed a method for inverting layer velocity using seismic wave propagation time (WT inversion). This method assumes that seismic waves follow the laws of sound wave propagation during their propagation in underground strata. Then, it uses the difference between the actual seismic wave travel time and the seismic wave travel time obtained from the calculated velocity model to invert the layer velocity. The velocity model is corrected multiple times through inversion calculation. When the difference between the actual seismic wave travel time and the calculated travel time meets the error requirement, it is considered that the inverted layer velocity is in good agreement with the actual layer velocity of the strata, thus achieving the purpose of layer velocity inversion. This method combines the advantages of conventional travel time inversion and full waveform inversion. It makes no high-frequency assumptions, and the algorithm converges quickly even if the initial velocity model differs significantly from the actual velocity. It also avoids assumptions similar to the Born and Rytov approximations. Furthermore, the algorithm remains unaffected by strong noise in the data and can invert the approximate outline of complex geological structures, thus possessing broad applicability. All these advantages make the WT layer velocity inversion method superior to other layer velocity inversion methods. However, the WT layer velocity inversion method calculates the gradient of the steepest descent direction using time difference backpropagation for all grid points in the velocity model, and velocity iteration is also performed on every grid point in the entire velocity model. This results in points not belonging to the transmitted wave propagation path being processed as well, leading to a significant increase in computational workload and a noticeable impact on inversion accuracy. Therefore, to improve the efficiency of the WT layer velocity inversion method, it is necessary to refine it. Summary of the Invention

[0004] In view of this, this application proposes a layer velocity inversion method based on the propagation law of sound waves in underground media, which only performs inversion calculations on grid points belonging to the sound wave propagation path. This method is the same in principle as the WT layer velocity inversion method, except that the range of the inversion grid is limited. Experimental results show that compared with the WT layer velocity inversion method, this method can greatly reduce the amount of computation and also greatly improve the accuracy of layer velocity inversion.

[0005] In a first aspect, embodiments of this application provide an improved method for inverting layer velocities using first-arrival waves, including:

[0006] Step S1: Obtain the result from X s Point excitation, X r The measured wave field value p(X) received at time t. r ,t;X s ) obs The measured wave field value p(X) was calculated. r ,t;X s ) obs derivative And pick up the measured wave field value p(X) r ,t;X s ) obs The first arrival time t(X) of the corresponding earthquake record r ,X s ) obs ;

[0007] Step S2: Determine the initial velocity model x(X)0;

[0008] Step S3: Obtain the velocity calculated using the current velocity model from X s Point excitation, X r The calculated wave field value p(X) received at time t. r ,t;X s ) cal And the calculated wave field value p(X) is obtained. r ,t;X s ) cal derivative

[0009] Step S4: Calculate the excitation point X s and receiving point X r The initial arrival time t(X,X) of each grid point in the velocity model s ) cal and t(X,X) r ) cal ;

[0010] Step S5: Based on the initial arrival time t(X,X) s ) cal and t(X,X) r ) cal To determine the propagation path of the first arrival wave;

[0011] Step S6: Calculate the new velocity model;

[0012] Step S7: Determine whether the new velocity model meets the accuracy requirements. If it does, the current velocity model is the inverted velocity model. If it does not meet the accuracy requirements, repeat steps S3-S6 based on the current velocity model until the new velocity model meets the accuracy requirements.

[0013] In one possible implementation, The calculation formula is as follows:

[0014]

[0015] In one possible implementation, The calculation formula is as follows:

[0016]

[0017] In one possible implementation, in step S4, the excitation point X is calculated using a functional equation. s and receiving point X r The initial arrival time t(X,X) of each grid point in the velocity model s ) cal and t(X,X) r ) cal .

[0018] In one possible implementation, in step S5, the propagation path of the first arrival wave is determined using the reciprocity principle and Fermat's principle.

[0019] In one possible implementation, step S5 specifically involves: adding up the first arrival times from the receiving point and the excitation point to each grid point on the same depth layer, and finding the grid point with the smallest value to be a point on the wave field propagation path on that depth layer; then combining the points on the wave field propagation path on each depth layer to obtain the propagation path of the first arrival wave.

[0020] In one possible implementation, step S6 specifically includes:

[0021] Step S6.1: Calculate Δτ(X) r ,X s The specific formula is:

[0022] Δτ(X r ,X s )=t(X,X s ) obs -t(X,X s ) cal

[0023] Where t(X,X) s ) obs The measured excitation point X sThe initial arrival time of each grid point in the velocity model, t(X,X) s ) cal The calculated excitation point X s The initial arrival time to each grid point in the velocity model;

[0024] Step S6.2: Calculate r(X) k The specific formula is as follows:

[0025]

[0026] in, For X s Point excitation, X r The derivative of the measured wave field value received at time t+τ at the point; To calculate the value of X using the velocity model s Point excitation, X r The derivative of the wave field value received at time t at the point; For X s Point excitation, X r The derivative of the measured wave field value received at time t+Δτ at the point; Let c(X) be the time derivative of the Green's function; c(X) is the formation velocity at X.

[0027] Step S6.3: According to c(X) k+1 =c(X) k +α k ·r(X) k α k Using the step size, a new velocity model is obtained.

[0028] Secondly, embodiments of this application provide an improved device for inverting layer velocities using first-arrival waves, comprising:

[0029] The data acquisition module is used to acquire data from X. s Point excitation, X r The measured wave field value p(X) received at time t. r ,t;X s ) obs The measured wave field value p(X) was calculated. r ,t;X s ) obs derivative And pick up the measured wave field value p(X) r ,t;X s ) obs The first arrival time t(X) of the corresponding earthquake record r ,X s ) obs ;

[0030] The initial layer velocity model construction module is used to determine the initial layer velocity model c(X)0;

[0031] The wave field value calculation module is used to obtain the wave field value calculated using the current velocity model from X. s Point excitation, X r The calculated wave field value p(X) received at time t. r ,t;X s ) cal And the calculated wave field value p(X) is obtained. r ,t;X s ) cal derivative

[0032] The first arrival time calculation module is used to calculate the excitation point X. s and receiving point X r The initial arrival time t(X,X) of each grid point in the velocity model s ) cal and t(X,X) r ) cal ;

[0033] The propagation path determination module is used to determine the propagation path based on the initial arrival time t(X,X). s ) cal and t(X,X) r ) cal To determine the propagation path of the first arrival wave;

[0034] A new velocity module calculation module is used to calculate the new velocity model;

[0035] The determination module is used to determine whether the new velocity model meets the accuracy requirement. If the accuracy requirement is met, the current velocity model is the inverted velocity model. If the accuracy requirement is not met, the new velocity model is repeatedly built iteratively based on the current velocity model until the new velocity model meets the accuracy requirement.

[0036] Thirdly, embodiments of this application provide an electronic device, including:

[0037] processor;

[0038] Memory;

[0039] And a computer program, wherein the computer program is stored in the memory, the computer program including instructions that, when executed by the processor, cause the electronic device to perform the method described in any one of the first aspects.

[0040] Fourthly, embodiments of this application provide a computer-readable storage medium including a stored program, wherein, when the program is executed, it controls the device where the computer-readable storage medium is located to perform the method described in any one of the first aspects.

[0041] The layer velocity inversion method proposed in this application is based on the WT layer velocity inversion theory and is an improvement upon it. It utilizes the first arrival wavefield for layer velocity inversion and is primarily applied in inter-well tomographic seismic exploration and non-zero bias VSP seismic exploration. Numerical simulation results show that, compared to the WT layer velocity inversion technique, the proposed method significantly improves the accuracy of the inverted velocity model, resulting in a substantial increase in velocity resolution. The thickness of the transition layer in the inverted velocity model is significantly less than that obtained using the WT method. Furthermore, the proposed method exhibits strong anti-interference capabilities, effectively suppressing noise interference with seismic signals. Even when the initial velocity model deviates significantly from the actual situation, the algorithm converges quickly and can accurately invert the morphology of complex geological structures. In summary, the proposed layer velocity inversion method absorbs all the advantages of the WT method, but its inverted layer velocity resolution is significantly improved compared to the results obtained using the WT method. Attached Figure Description

[0042] To more clearly illustrate the technical solutions of the embodiments of this application, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0043] Figure 1 A schematic flowchart illustrating the improved method for inverting layer velocities using first-arrival waves, provided in an embodiment of this application;

[0044] Figure 2 A schematic diagram showing the high-speed anomaly horizontal interlayer velocity model provided in this application embodiment and the results of different inversion methods using the initial velocity as the background velocity;

[0045] Figure 3 A schematic diagram showing the fault layer velocity model perpendicular to the fault plane provided in the embodiments of this application, and the results of different inversion methods using the initial velocity as the background velocity;

[0046] Figure 4 A structural block diagram of an improved first-arrival inversion layer velocity device provided in this application embodiment;

[0047] Figure 5This is a schematic diagram of the structure of an electronic device provided in an embodiment of this application. Detailed Implementation

[0048] To better understand the technical solution of this application, the embodiments of this application will be described in detail below with reference to the accompanying drawings.

[0049] It should be understood that the described embodiments are merely some, not all, of the embodiments in this application. All other embodiments obtained by those skilled in the art based on the embodiments in this application without inventive effort are within the scope of protection of this application.

[0050] The terminology used in the embodiments of this application is for the purpose of describing particular embodiments only and is not intended to be limiting of this application. The singular forms “a,” “the,” and “the” used in the embodiments of this application and the appended claims are also intended to include the plural forms unless the context clearly indicates otherwise.

[0051] It should be understood that the term "and / or" used in this article is merely a description of the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent: A existing alone, A and B existing simultaneously, or B existing alone. Additionally, the character " / " in this article generally indicates that the preceding and following related objects have an "or" relationship.

[0052] Formation layer velocity is a crucial physical property parameter in seismic exploration technology. Obtaining accurate subsurface layer velocities has been a subject of ongoing exploration and research in seismic data processing and interpretation. Due to the widespread application of surface seismic technology in oil and gas exploration over the past few decades, researchers typically extracted subsurface layer velocities using first-order reflection wave information. In recent years, with the increasing application of non-zero biased VSP exploration technology and inter-well tomography, the goal of researchers has become how to invert subsurface layer velocities using seismic transmission wave fields. This invention improves upon the previously proposed WT layer velocity inversion method by utilizing the propagation path of seismic waves. Specifically, it performs layer velocity inversion only on velocity grid points within the seismic wave propagation path in the velocity model, replacing the WT method's method of inverting layer velocities for all grid points in the entire velocity model. The calculation results of the numerical velocity model show that, compared with the WT layer velocity inversion method, the layer velocity inversion method used in this invention can not only greatly reduce the amount of calculation and shorten the time for layer velocity inversion, but also greatly improve the accuracy of layer velocity inversion, making the inverted layer velocity model closer to the actual layer velocity model.

[0053] The layer velocity inversion method based on transmitted wave fields proposed in this application is an improvement on the WT layer velocity inversion method by utilizing the propagation path of transmitted wave fields. Its layer velocity inversion principle is the same as the WT layer velocity inversion method. However, unlike the WT method, the layer velocity inversion method used in this application constrains the scope of layer velocity inversion; that is, it only inverts velocity grid points belonging to the wave field propagation path. Therefore, it is a constrained WT layer velocity inversion. Thus, the layer velocity inversion method used in this invention can be divided into the following two parts:

[0054] (I) WT layer velocity inversion

[0055] Assume that the propagation of seismic waves satisfies the following two-dimensional acoustic wave equation:

[0056]

[0057] Where c(X) represents the formation velocity at X, and p(X,t; X s ) represents the pressure wave field value at X, ρ(X) represents the formation density at X, and S(t;X) represents the pressure wave field value at X. s ) represents the source function.

[0058] p(X) r ,t;X s ) obs X represents s Point excitation, X r The measured wavefield value received at time t by the point, p(X) r ,t;X s ) cal X represents s Point excitation, X r The wave field value received at time t using the velocity model, and the mediation function f(X) r ,τ;X s Let ) represent the cross-correlation function between the measured wavefield and the calculated wavefield, then we have:

[0059] f(X r ,τ;X s )=∫p(X r ,t+τ;X s ) obs p(X r ,t;X s ) cal dt

[0060] Where p(X) r ,t+τ;X s ) obs X represents s Point excitation, X rThe measured wave field value received at time t+τ, where τ represents the translation factor.

[0061] To ensure that the measured wave field value p(X) r ,t;X s ) obs The calculated wave field value p(X) r ,t;X s ) cal To achieve the best match, we need to calculate an optimal translation factor τ. In other words, we need to find a Δτ such that f(X) r ,τ;X s The value of f(X) is the largest, and according to the extremum property of a function, f(X) is the largest. r ,τ;X s The derivative with respect to τ should be 0 when τ = Δτ, that is, the following formula holds:

[0062]

[0063] in,

[0064]

[0065] Determine the mediation function f(X) r ,τ;X s After that, the mismatch function S is defined as follows:

[0066]

[0067] Where s and r are the excitation point number and the receiving point number, respectively, and X r Indicates the receiving point, X s Indicates the excitation point.

[0068] In the velocity inversion of the WT layer, the gradient of the steepest descent direction of the mismatch function S is used as a measure of changing the computational velocity model, assuming c(X). k Given the velocity model after k iterations, the mismatch function S with respect to velocity c(X) is obtained. k The derivative r(X) k , r(X) k The expression is as follows:

[0069]

[0070] Then, the velocity model c(X) obtained after the (k+1)th iteration is... k+1 as follows:

[0071] c(X) k+1 =c(X) k +α k ·r(X) k

[0072] Where, α k The step size.

[0073] Through derivation, we obtain r(X). k The final expression is as follows:

[0074]

[0075] in,

[0076]

[0077] Let be the time derivative of the Green's function.

[0078] (II) Determination of the wave field propagation path

[0079] In the stratigraphic velocity inversion method proposed in this application, determining the seismic wave propagation path requires the application of the acoustic reciprocity principle and Fermat's principle. The acoustic reciprocity principle states that the propagation path of a seismic wave from the source to the receiver is the same as the path from the receiver to the source. Fermat's principle states that waves propagate along the path with the shortest travel time between two points. Therefore, by summing the first arrival times of all velocity grid points at the same depth from the receiver and excitation points, the grid point with the smallest sum belongs to the wave field propagation path. The combination of such grid points at each depth layer forms the propagation path of the wave field. In the stratigraphic velocity inversion method proposed in this application, the first arrival times of each grid point in the velocity model are calculated using the finite difference approximation equation. Using the equation to calculate the first arrival times has the advantages of no blind spots, low computational cost, and high accuracy. In a two-dimensional isotropic medium, the mathematical expression of the equation is as follows:

[0080]

[0081] Where x and z are the horizontal and depth coordinates of the grid point, respectively, and s(x,z) is the slowness distribution function of the velocity model (i.e., the reciprocal of the velocity).

[0082] The following is a detailed explanation with reference to the accompanying drawings.

[0083] See Figure 1 This is a schematic flowchart illustrating the improved method for retrieving layer velocities using first-arrival waves, provided in an embodiment of this application. Figure 1 As shown, it mainly includes the following steps.

[0084] Step S1: Obtain the result from X s Point excitation, X r The measured wave field value p(X) received at time t. r,t;X s ) obs The measured wave field value p(X) was calculated. r ,t;X s ) obs derivative And pick up the measured wave field value p(X) r ,t;X s ) obs The first arrival time t(X) of the corresponding earthquake record r ,X s ) obs .

[0085] Step S2: Determine the initial velocity model c(X)0.

[0086] Step S3: Obtain the velocity calculated using the current velocity model from X s Point excitation, X r The calculated wave field value p(X) received at time t. r ,t;X s ) cal And the calculated wave field value p(X) is obtained. r ,t;X s ) cal derivative

[0087] Step S4: Use the equation to find the excitation point X. s and receiving point X r The initial arrival time t(X,X) of each grid point in the velocity model s ) cal and t(X,X) r ) cal ;

[0088] Step S5: Based on the initial arrival time t(X,X) s ) cal and t(X,X) r ) cal The propagation path of the first arrival wave is determined using the reciprocity principle and Fermat's principle. Specifically:

[0089] The first arrival times from the receiving point and the excitation point to each grid point at the same depth are added together. The grid point with the smallest value belongs to the point on the wave field propagation path at that depth. The propagation path of the first arrival wave is obtained by combining the points on the wave field propagation path at each depth.

[0090] Step S6: Calculate the new velocity model. Specifically:

[0091] Step S6.1: Calculate Δτ(X) r ,X sThe specific formula is:

[0092] Δτ(X r ,X s )=t(X,X s ) obs -t(X,X s ) cal

[0093] Step S6.2: Calculate r(X) k The specific formula is as follows:

[0094]

[0095] in, Let be the time derivative of the Green's function.

[0096] Step S6.3: According to c(X) k+1 =c(X) k +α k ·r(X) k α k Using the step size, a new velocity model is obtained.

[0097] Step S7: Determine whether the new velocity model meets the accuracy requirements. If it does, the current velocity model is the inverted velocity model. If it does not meet the accuracy requirements, repeat steps S3-S6 based on the current velocity model until the new velocity model meets the accuracy requirements.

[0098] Step S7: Determine whether the new velocity model meets the accuracy requirements. If it does, the current velocity model is the inverted velocity model. If it does not meet the accuracy requirements, repeat steps S3-S6 based on the current velocity model until the new velocity model meets the accuracy requirements.

[0099] According to the technical solution provided in the embodiments of this application, two numerical layer velocity models are selected for numerical simulation calculations. A constant velocity model with an initial velocity equal to the background velocity is selected as the initial velocity model. The selected excitation observation system is an inter-well tomographic seismic observation system. The excitation sources and receivers are uniformly distributed to the right of the velocity model. It is assumed that the grid is divided into 2.4m*2.4m sections, the time sampling interval is 0.25ms, and the source function used is a Ricker wavelet with a main frequency of 80Hz. First, the forward wavefield records of the actual velocity model and the initial velocity model are calculated using the finite difference method of the acoustic wave equation. Then, inversion is performed using the WT layer velocity inversion method and the layer velocity inversion method proposed in this invention, respectively. After multiple iterations and calculations, the inverted layer velocity models are finally obtained.

[0100] See Figure 2This is a schematic diagram illustrating the high-speed anomaly horizontal interlayer velocity model provided in this application embodiment, and the results of different inversion methods using the initial velocity as the background velocity. Figure 2 As shown, Figure 2 (a) shows a layer velocity model with a horizontal high-speed anomaly interlayer. The background velocity is 2000 m / s, the layer velocity of the horizontal interlayer is 2500 m / s, and the thickness of the horizontal interlayer is 14.4 m. Figure 2 (b) is the result of 12 iterations of the WT layer velocity inversion method; Figure 2 (c) The result of 12 iterations of the layer velocity inversion method proposed in this application. Comparing the inversion results of the two inversion methods, it is easy to see that both the WT layer velocity inversion and the inversion method proposed in this application can correctly invert the location of the high-speed anomaly interlayer. However, the layer velocity inversion method of this application can more accurately reflect the velocity characteristics of the anomaly, and the inversion result is closer to the actual velocity model.

[0101] See Figure 3 This is a schematic diagram illustrating the fault layer velocity model perpendicular to the fault plane provided in this application embodiment, and the results of different inversion methods using the initial velocity as the background velocity. Figure 3 As shown, Figure 3 (a) shows a fault velocity model perpendicular to the fault plane, with a background velocity of 2000 m / s and a fault anomaly velocity of 2500 m / s. Figure 3 (b) is the result of 12 iterations of the WT layer velocity inversion method; Figure 3 (c) The result of 12 iterations of the layer velocity inversion method proposed in this application. Comparing the inversion results of the two inversion methods, it is easy to see that both the WT layer velocity inversion and the inversion method proposed in this application can accurately invert the location of the uplift and downlift blocks of the fault. However, the layer velocity inversion method of this application is significantly more accurate than the WT layer velocity inversion method.

[0102] In summary, the layer velocity inversion method proposed in this application is based on the WT layer velocity inversion theory and is an improvement upon it. It utilizes the first arrival wave field for layer velocity inversion and is primarily applied in inter-well tomographic seismic exploration and non-zero bias VSP seismic exploration. Numerical simulation results show that, compared to the WT layer velocity inversion technique, the proposed method significantly improves the accuracy of the inverted velocity model, resulting in a substantial increase in velocity resolution. The thickness of the transition layer in the inverted velocity model is significantly smaller than that obtained using the WT method. Furthermore, the proposed method exhibits strong anti-interference capabilities, effectively suppressing noise interference with seismic signals. Even when the initial velocity model deviates significantly from the actual situation, the algorithm converges quickly, accurately reversing the morphology of complex geological structures. In general, the proposed method absorbs all the advantages of the WT method, but its inverted layer velocity resolution is significantly improved compared to the results obtained using the WT method.

[0103] Corresponding to the above embodiments, this application also provides an improved device for inverting layer velocities using first arrival waves.

[0104] See Figure 4 This is a structural block diagram of an improved first-arrival wave inversion layer velocity device provided in an embodiment of this application. Figure 4 As shown, it mainly includes the following modules.

[0105] Data acquisition module 401 is used to acquire data from X s Point excitation, X r The measured wave field value p(X) received at time t. r ,t;X s ) obs The measured wave field value p(X) was calculated. r ,t;X s ) obs derivative And pick up the measured wave field value p(X) r ,t;X s ) obs The first arrival time t(X) of the corresponding earthquake record r ,X s ) obs .

[0106] The initial layer velocity model construction module 402 is used to determine the initial layer velocity model c(X)0.

[0107] Wave field value calculation module 403 is used to obtain the wave field value calculated using the current velocity model from X. s Point excitation, Xr The calculated wave field value p(X) received at time t. r ,t;X s ) cal And the calculated wave field value p(X) is obtained. r ,t;X s ) cal derivative

[0108] First arrival time calculation module 404 is used to calculate the excitation point X using the equation of process function. s and receiving point X r The initial arrival time t(X,X) of each grid point in the velocity model s ) cal and t(X,X) r ) cal .

[0109] The propagation path determination module 405 is used to determine the propagation path based on the initial arrival time t(X,X). s ) cal and t(X,X) r ) cal The propagation path of the first arrival wave is determined by using the reciprocity principle and Fermat's principle.

[0110] The new velocity module calculation module 406 is used to calculate the new velocity model.

[0111] The determination module 407 is used to determine whether the new velocity model meets the accuracy requirement. If the accuracy requirement is met, the current velocity model is the inverted velocity model. If the accuracy requirement is not met, the new velocity module is repeatedly built iteratively based on the current velocity model until the new velocity model meets the accuracy requirement.

[0112] It should be noted that the specific content involved in the embodiments of this application can be found in the description of the above method embodiments, and will not be repeated here for the sake of brevity.

[0113] Corresponding to the above embodiments, this application also provides an electronic device.

[0114] See Figure 5 This is a schematic diagram of the structure of an electronic device provided in an embodiment of this application. Figure 5 As shown, the electronic device 500 may include a processor 501, a memory 502, and a communication unit 503. These components communicate via one or more buses. Those skilled in the art will understand that the electronic device structure shown in the figure does not constitute a limitation on the embodiments of this application. It may be a bus topology or a star topology, and may include more or fewer components than shown, or combine certain components, or have different component arrangements.

[0115] The communication unit 503 is used to establish a communication channel, thereby enabling the electronic device to communicate with other devices.

[0116] Processor 501 serves as the control center of the electronic device, connecting various parts of the device via various interfaces and lines. It executes software programs and / or modules stored in memory 502, and calls data stored in memory to perform various functions and / or process data. The processor may be composed of integrated circuits (ICs), such as a single packaged IC or multiple packaged ICs with the same or different functions connected together. For example, processor 501 may consist only of a central processing unit (CPU). In this embodiment, the CPU may have a single processing core or include multiple processing cores.

[0117] Memory 502 is used to store the execution instructions of processor 501. Memory 502 can be implemented by any type of volatile or non-volatile storage device or a combination thereof, such as static random access memory (SRAM), electrically erasable programmable read-only memory (EEPROM), erasable programmable read-only memory (EPROM), programmable read-only memory (PROM), read-only memory (ROM), magnetic storage, flash memory, magnetic disk or optical disk.

[0118] When the execution instructions in memory 502 are executed by processor 501, the electronic device 500 is able to perform some or all of the steps in the above method embodiments.

[0119] Corresponding to the above embodiments, this application also provides a computer-readable storage medium, wherein the computer-readable storage medium may store a program, wherein when the program runs, it can control the device where the computer-readable storage medium is located to execute some or all of the steps in the above method embodiments. Specifically, the computer-readable storage medium may be a magnetic disk, an optical disk, read-only memory (ROM), or random access memory (RAM), etc.

[0120] Corresponding to the above embodiments, this application also provides a computer program product containing executable instructions that, when executed on a computer, cause the computer to perform some or all of the steps in the above method embodiments.

[0121] In this application embodiment, "at least one" refers to one or more, and "more than one" refers to two or more. "And / or" describes the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent the existence of A alone, the simultaneous existence of A and B, or the existence of B alone. A and B can be singular or plural. The character " / " generally indicates that the preceding and following related objects are in an "or" relationship. "At least one of the following" and similar expressions refer to any combination of these items, including any combination of single or plural items. For example, at least one of a, b, and c can represent: a, b, c, ab, ac, bc, or abc, where a, b, and c can be single or multiple.

[0122] Those skilled in the art will recognize that the units and algorithm steps described in the embodiments disclosed herein can be implemented using electronic hardware, computer software, or a combination of electronic hardware and software. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0123] Those skilled in the art will understand that, for the sake of convenience and brevity, the specific working processes of the systems, devices, and units described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0124] In the several embodiments provided in this application, any function, if implemented as a software functional unit and sold or used as an independent product, can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, or a part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of this application. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0125] The above description is merely a specific embodiment of this application. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the protection scope of this application. The protection scope of this application should be determined by the protection scope of the claims.

Claims

1. An improved method for inverting layer velocities using first-arrival waves, characterized in that, include: Step S1: Obtain the result from X s Point excitation, X r The measured wave field value p(X) received at time t. r ,t;X s ) obs The measured wave field value p(X) was calculated. r ,t;X s ) obs derivative And pick up the measured wave field value p(X) r ,t;X s ) obs The first arrival time t(X) of the corresponding earthquake record r ,X s ) obs ; Step S2: Determine the initial velocity model c(X)0; Step S3: Obtain the velocity calculated using the current velocity model from X s Point excitation, X r The calculated wave field value p(X) received at time t. r ,t;X s ) cal And the calculated wave field value p(X) is obtained. r ,t;X s ) cal derivative Step S4: Calculate the excitation point X s and receiving point X r The initial arrival time t(X,X) of each grid point in the velocity model s ) cal and t(X,X) r ) cal ; Step S5: Based on the initial arrival time t(X,X) s ) cal and t(X,X) r ) cal To determine the propagation path of the first arrival wave; Step S6: Calculate the new velocity model; Step S7: Determine whether the new velocity model meets the accuracy requirements. If it does, the current velocity model is the inverted velocity model. If it does not meet the accuracy requirements, repeat steps S3-S6 based on the current velocity model until the new velocity model meets the accuracy requirements.

2. The improved method for inverting layer velocities using first-arrival waves according to claim 1, characterized in that, The calculation formula is as follows:

3. The improved method for inverting layer velocities using first-arrival waves according to claim 1, characterized in that, The calculation formula is as follows:

4. The improved method for inverting layer velocities using first-arrival waves according to claim 1, characterized in that, In step S4, the excitation point X is calculated using the equation of process function. s and receiving point X r The initial arrival time t(X,X) of each grid point in the velocity model s ) cal and t(X,X) r ) cal .

5. The improved method for inverting layer velocities using first-arrival waves according to claim 1, characterized in that, In step S5, the propagation path of the first arrival wave is determined using the reciprocity principle and Fermat's principle.

6. The improved method for inverting layer velocities using first-arrival waves according to claim 1, characterized in that, Step S5 specifically involves: adding up the first arrival times from the receiving point and the excitation point to each grid point at the same depth of the stratum, and finding the grid point with the smallest value to be a point on the wave field propagation path at that depth of the stratum. Then, combining the points on the wave field propagation path at each depth of the stratum, the propagation path of the first arrival wave is obtained.

7. The improved method for inverting layer velocities using first-arrival waves according to claim 1, characterized in that, Step S6 specifically involves: Step S6.1: Calculate Δτ(X) r ,X s The specific formula is: Δτ(X r ,X s )=t(X,X s ) obs -t(X,X s ) cal Where t(X,X) s ) obs The measured excitation point X s The initial arrival time of each grid point in the velocity model, t(X,X) s ) cal The calculated excitation point X s The initial arrival time to each grid point in the velocity model; Step S6.2: Calculate r(X) k The specific formula is as follows: in, For X s Point excitation, X r The derivative of the measured wave field value received at time t+τ at the point; To calculate the value of X using the velocity model s Point excitation, X r The derivative of the wave field value received at time t at the point; For X s Point excitation, X r The derivative of the measured wave field value received at time t+Δτ at the point; Let c(X) be the time derivative of the Green's function; c(X) is the formation velocity at X. Step S6.3: According to c(X) k+1 =c(X) k +α k ·r(X) k α k Using the step size, a new velocity model is obtained.

8. An improved device for inverting layer velocity using first-arrival waves, characterized in that, include: The data acquisition module is used to acquire data from X. s Point excitation, X r The measured wave field value p(X) received at time t. r ,t;X s ) obs The measured wave field value p(X) was calculated. r ,t;X s ) obs derivative And pick up the measured wave field value p(X) r ,t;X s ) obs The first arrival time t(X) of the corresponding earthquake record r ,X s ) obs ; The initial layer velocity model construction module is used to determine the initial layer velocity model c(X)0; The wave field value calculation module is used to obtain the wave field value calculated using the current velocity model from X. s Point excitation, X r The calculated wave field value p(X) received at time t. r ,t;X s ) cal And the calculated wave field value p(X) is obtained. r ,t;X s ) cal derivative The first arrival time calculation module is used to calculate the excitation point X. s and receiving point X r The initial arrival time t(X,X) of each grid point in the velocity model s ) cal and t(X,X) r ) cal ; The propagation path determination module is used to determine the propagation path based on the initial arrival time t(X,X). s ) cal and t(X,X) r ) cal To determine the propagation path of the first arrival wave; A new velocity module calculation module is used to calculate the new velocity model; The determination module is used to determine whether the new velocity model meets the accuracy requirement. If the accuracy requirement is met, the current velocity model is the inverted velocity model. If the accuracy requirement is not met, the new velocity model is repeatedly built iteratively based on the current velocity model until the new velocity model meets the accuracy requirement.

9. An electronic device, characterized in that, include: processor; Memory; And a computer program, wherein the computer program is stored in the memory, the computer program including instructions that, when executed by the processor, cause the electronic device to perform the method of any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium includes a stored program, wherein, when the program is executed, it controls the device on which the computer-readable storage medium is located to perform the method according to any one of claims 1 to 7.

Images