Fast beam shift imaging method, device and apparatus

Abstract

The present application provides a fast beam migration imaging method, apparatus, and device, comprising: dividing seismic data of a target area into a series of supergathers; for each supergather, identifying the arrival time of a local event based on its energy; determining the slope direction of the local event at the arrival time position; superimposing the local events into wavelets along the slope direction of the local events to construct a beam; performing kinematic and dynamic ray tracing along the slope direction of the corresponding local event at the center position of the supergather to determine the imaging position of the beam; at the imaging position of the beam, expanding the beam based on the real travel time and imaginary travel time of a Gaussian beam from the center position of the supergather; and accumulating the beam expansion results of each supergather to output an imaging result. The present application can improve the imaging effect of complex structural areas while addressing the calculation problem of cross-event slope.

Description

Technical Field

[0001] The present application relates to the field of seismic migration imaging technology, and in particular to a rapid beam migration imaging method, apparatus, and device. Background Art

[0002] The rapid beam migration imaging method decomposes seismic records into a series of beams (wavelets with specific arrival times, shot and receiver locations, and shot and receiver slope directions). Imaging then requires only the defined beam directions (shot and receiver slope directions). Beam decomposition only needs to be performed once (independent of velocity), resulting in high speed and precision. Rapid beam migration primarily involves three steps: beam decomposition, beam positioning, and beam expansion. Accurate beam decomposition is a prerequisite for accurate beam positioning and expansion, and plays a decisive role in the final imaging quality. The most critical step in beam decomposition is determining the local event slope direction (the shot and receiver slope directions). The accuracy of this local event slope direction directly impacts the accuracy of beam decomposition.

[0003] In the prior art, there are several methods for obtaining the slope direction of the local event axis: the first is to use the local tilt stacking method to obtain it in multiple conventional gathers (common shot, common receiving point and common shot offset gathers) at the same time; the second is to use the multidimensional local tilt stacking method to obtain it directly in the super gather; and the third is to use the plane wave decomposition filter method to obtain it in the super gather.

[0004] None of the aforementioned techniques effectively handles events crossing. However, this phenomenon is very common in seismic data, especially in complex tectonic regions (such as unconformities or faults). This problem hinders fast beam migration imaging in complex tectonic regions. Furthermore, determining slopes in supergathers requires simultaneous determination of slopes in four directions in a three-dimensional context, which creates severe high-dimensional computational challenges and compromises the efficiency of beam decomposition. Summary of the Invention

[0005] In view of the problems in the prior art that the situation of crossed phase axes cannot be handled well and the beam decomposition calculation efficiency is low, resulting in poor imaging effects in complex structural areas, the present application aims to provide a fast beam offset imaging method, device and equipment, which can improve the imaging efficiency and imaging quality in complex structural areas.

[0006] In a first aspect, the present application provides a fast beam shift imaging method, comprising:

[0007] Divide the seismic data of the target area into a series of super gathers;

[0008] For each super gather, the arrival time of the local event is identified based on the energy of the local event at the center of the super gather;

[0009] Determining the slope direction of the local event at the arrival time position of the local event, wherein the slope direction includes the shot point slope direction and the detection point slope direction;

[0010] The local events are superimposed into wavelets along the slope direction of the local events to construct a beam;

[0011] Perform kinematic and dynamic ray tracing along the slope direction of the corresponding local event at the center of the beam super gather, and determine the imaging position of the beam according to the temporal and spatial positioning conditions;

[0012] At the imaging position of the beam, the beam is expanded according to the real travel time and virtual travel time of the Gaussian beam at the center of the super gather to obtain the beam expansion result;

[0013] The beam expansion results of each super gather are accumulated to output the imaging result.

[0014] In a possible embodiment, determining the slope direction of the local event according to the arrival time of the local event includes:

[0015] Determining the slope direction of the local event at the arrival time position of the local event includes:

[0016] The differential evolution algorithm based on neighborhood species is used to drive the distribution of individuals in the population to the vicinity of the local optimal solution;

[0017] The population individuals near the local optimal solution are screened by the distance-based individual screening algorithm to obtain the excellent individuals closest to each local optimal solution;

[0018] The Hooke-Jeeves pattern search algorithm is used to ensure that the selected excellent individuals converge to the local optimal solution with certainty.

[0019] In a possible embodiment, it further includes:

[0020] The obtained local optimal solutions are further screened using a distance-based individual screening algorithm.

[0021] In a possible embodiment, dividing the seismic data of the target area into a series of super gathers includes:

[0022] Shot lines are defined according to the start and end line numbers and line intervals of the vertical and horizontal survey lines of the observation system. The intersection point on any shot line is used as the center position of the shot point of the super gather. All shot points within the first preset range are selected as the super gather.

[0023] The detection point lines are defined according to the vector offset range and the offset sampling interval, and the intersection of any detection point line is used as the detection point center position of the super gather. All detection points within the second preset range are selected as the super gather.

[0024] In a possible embodiment, for each super gather, identifying the arrival time of the local event according to the energy of the local event at the center of the super gather includes:

[0025] The method of identifying the arrival time of the local event according to the energy of the local event at the center of each super gather comprises:

[0026] Calculate the energy curve of the local event in a single super gather;

[0027] determining a local maximum value of an energy curve according to the energy curve of the local event;

[0028] The arrival time of the local event is determined according to the local maximum value of the energy curve.

[0029] In a possible embodiment, performing kinematic and dynamic ray tracing along the slope direction of the corresponding local event at the center position of the super gather of the beam and determining the imaging position of the beam according to time and space positioning conditions includes:

[0030] Solving the kinematic ray tracing equation and the dynamic ray tracing equation to obtain the ray paths from the shot point and the receiver point respectively and the complex travel time of the Gaussian beam around the ray path;

[0031] The point closest to the arrival time of the local event is found near the closest position between the ray paths from the shot point and the receiver point, and is used as the imaging position of the beam.

[0032] In a possible embodiment, the step of expanding the beam at the imaging position of the beam according to the real travel time and the virtual travel time of the Gaussian beam at the center position of the super gather to obtain the beam expansion result includes:

[0033] For any imaging point near the imaging position of the beam, determine the sum of the real part and the sum of the imaginary part of the complex travel time of the Gaussian beam from the shot point and the receiver point;

[0034] The image field value of the imaging point is obtained by interpolating the sum of the real parts of the Gaussian beam in the wavelet of the beam;

[0035] The image field value is attenuated by utilizing the sum of the imaginary parts of the Gaussian beam, thereby obtaining a beam expansion result.

[0036] In a second aspect, the present application further provides a rapid beam shift imaging device, comprising:

[0037] A data processing module for dividing the seismic data of the target area into a series of super gathers;

[0038] An arrival time identification module is used to identify the arrival time of the local event based on the energy of the local event at the center of each super gather;

[0039] A slope direction determination module is used to determine the slope direction of the local event at the arrival time position of the local event, wherein the slope direction includes the shot point slope direction and the detection point slope direction;

[0040] A beam decomposition module is used to superimpose local events into wavelets along the slope direction of the local events to construct a beam;

[0041] The beam positioning module is used to perform kinematic and dynamic ray tracing along the slope direction of the corresponding local phase axis at the center position of the beam super gather, and determine the imaging position of the beam according to the time and space positioning conditions;

[0042] A beam spreading module is used to spread the beam at the imaging position of the beam according to the real travel time and the virtual travel time of the Gaussian beam at the center position of the super gather to obtain a beam spreading result;

[0043] The imaging output module is used to accumulate the beam expansion results of each super gather and output the imaging results.

[0044] In a third aspect, the present application further provides a rapid beam shift imaging system, comprising: at least one processor, and a memory communicatively connected to the at least one processor;

[0045] The memory stores computer-executable instructions;

[0046] The processor executes the computer-executable instructions stored in the memory to implement the method in any possible implementation of the first aspect above.

[0047] In a fourth aspect, the present application further provides a computer-readable storage medium, wherein the computer-readable storage medium stores computer-executable instructions, and when the computer-executable instructions are executed by a processor, they are used to implement the method in any possible implementation manner of the first aspect.

[0048] The fast beam migration imaging method, device and equipment provided by the present application are multi-peaked because the coherence surface of the cross-event axes along different slope directions in the super gather is multi-peaked, and the position of the local extreme value accurately corresponds to the slope direction of the corresponding event axis. Therefore, the present application converts the calculation problem of the slope of the cross-event axes into the problem of searching for multiple local extreme values ​​of the coherence surface. That is, the coherence of the event axes along different slope directions can be used as the objective function and converted into a multi-peak optimization problem. The present application solves the slope of the local event axis in the fast beam migration beam decomposition link through a differential evolution algorithm based on neighborhood species, an individual screening algorithm based on distance and a Hooke-Jeeves pattern search algorithm, thereby obtaining high-precision calculation results. This method can reduce the impact of high-dimensional calculation problems on computational efficiency in the case of three-dimensional super gathers, thereby improving the imaging effect of complex structural areas when dealing with the calculation problem of the cross-event axis slope. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] The accompanying drawings, which are incorporated in and constitute a part of this specification, illustrate embodiments consistent with the present application and, together with the description, serve to explain the principles of the present application.

[0050] Figure 1 A flowchart of a fast beam shift imaging method provided in one embodiment of the present application;

[0051] Figure 2 Schematic diagram of 2D and 3D super gather division;

[0052] Figure 3a This is a schematic diagram of 2D supergather cross-event seismic data;

[0053] Figure 3b for Figure 3a The coherence result diagram obtained by using the variance-based coherence calculation formula at the intersection of the event along different slope directions at the super channel concentration;

[0054] Figure 4 Schematic diagram for identifying the arrival time of local events in a two-dimensional super gather;

[0055] Figure 5 Schematic diagram for Hooke-Jeeves pattern search;

[0056] Figure 6 Schematic diagram of beam positioning;

[0057] Figure 7 Schematic diagram of beam expansion;

[0058] Figure 8 It is the Marmousi2 speed model;

[0059] Figure 9This is the migration result diagram based on the Marmousi2 velocity model;

[0060] Figure 10 This is a partial enlarged view of the Marmousi2 velocity model migration result;

[0061] Figure 11 Velocity profiles at different locations of the 3D SEG / EAGE salt model;

[0062] Figure 12 This is the migration result of the 3D SEG / EAGE salt model at z = 1020m;

[0063] Figure 13 The migration result of the 3D SEG / EAGE salt model at x = 4620m is shown;

[0064] Figure 14 The migration result diagram of the 3D SEG / EAGE salt model at x = 6260m;

[0065] Figure 15 This is the migration result of the 3D SEG / EAGE salt model at y = 6280m;

[0066] Figure 16 A schematic structural diagram of a fast beam shift imaging device provided in one embodiment of the present application;

[0067] Figure 17 A schematic structural diagram of a rapid beam shift imaging system provided in one embodiment of the present application.

[0068] The above drawings illustrate specific embodiments of the present application, which will be described in more detail below. These drawings and the textual description are not intended to limit the scope of the present application in any way, but rather to illustrate the concepts of the present application to those skilled in the art by reference to specific embodiments. DETAILED DESCRIPTION

[0069] Exemplary embodiments will be described in detail herein, with examples illustrated in the accompanying drawings. In the following description, when referring to the drawings, identical numerals in different figures represent identical or similar elements, unless otherwise indicated. The embodiments described in the following exemplary embodiments are not intended to represent all embodiments consistent with the present application. Rather, they are merely examples of apparatus and methods consistent with certain aspects of the present application, as detailed in the appended claims.

[0070] In the embodiments of the present application, words such as "first" and "second" are used to distinguish identical or similar items with substantially the same functions and effects. Those skilled in the art will understand that words such as "first" and "second" do not limit the quantity or execution order, and words such as "first" and "second" do not necessarily mean different.

[0071] It should be noted that, in the embodiments of the present application, words such as "exemplary" or "for example" are used to indicate examples, illustrations or descriptions. Any embodiment or design described in this application as "exemplary" or "for example" should not be interpreted as being more preferred or advantageous than other embodiments or design. Specifically, the use of words such as "exemplary" or "for example" is intended to present related concepts in a concrete way. In the embodiments of the present application, "at least one" refers to one or more, and "a plurality" refers to two or more.

[0072] It should be noted that the “at…” in the embodiments of the present application can be the instant when a certain situation occurs, or it can be a period of time after the occurrence of a certain situation. The embodiments of the present application do not make specific limitations on this.

[0073] Fast beam migration uses beam decomposition to decompose seismic records into a series of beams (wavelets with specific arrival times, shot and receiver locations, and shot and receiver slope directions). Imaging then requires only the defined beam directions (shot and receiver slope directions). Beam decomposition only needs to be performed once (independent of velocity), resulting in high speed and accuracy. Fast beam migration is computationally two orders of magnitude faster than Kirchhoff depth migration, commonly used in industrial production.

[0074] Rapid beam migration primarily involves three steps: beam decomposition, beam positioning, and beam expansion. Accurate beam decomposition is a prerequisite for accurate beam positioning and expansion, and plays a decisive role in the final imaging quality. The most critical step in beam decomposition is determining the local event slope direction (the slope direction of the shot and receiver points). The accuracy of this local event slope direction directly affects the accuracy of beam decomposition.

[0075] In the prior art, there are three main ways to obtain the local event slope direction in the fast beam migration beam decomposition link: the first is to use the local tilt stacking method to obtain it in multiple conventional gathers (common shot, common receiving point and common shot offset gathers) at the same time; the second is to use the multidimensional local tilt stacking method to obtain it directly in the super gather; and the third is to use the plane wave decomposition filter method to obtain it in the super gather.

[0076] However, the first method is affected by the large sampling interval and event crossing phenomenon along the transverse line in three-dimensional conditions, which is not only time-consuming but also difficult. The second method requires sampling along four directions (p sx 、p sy 、p rx 、p ry ) uniform sampling will cause serious high-dimensional calculation problems. For example, if the number of sampling points along each direction is set to 100, 10 8 The third method also needs to face the problem of high-dimensional calculation in three-dimensional cases, and the plane wave decomposition filter method is difficult to handle the phenomenon of event crossing. It can be seen that the method of obtaining the slope direction of the local event in the prior art cannot handle the situation of event crossing well. However, for seismic data, the phenomenon of event crossing is very common, especially for seismic data bodies from complex structural areas (such as unconformities or faults). Affected by this problem, the imaging effect of fast beam migration in complex structural areas is poor. In addition, if the slope is obtained in the super gather, the slopes in four directions need to be obtained simultaneously in the three-dimensional case, which will cause serious high-dimensional calculation problems, thereby affecting the computational efficiency of beam decomposition.

[0077] In order to solve the above problems, the present application provides a fast beam migration imaging method, device and equipment. Through the fast beam migration imaging method, the method can effectively reduce the impact of high-dimensional calculation problems on computational efficiency in the case of three-dimensional super gathers, and improve the imaging effect of complex structural areas when dealing with the calculation problem of cross-phase axis slope.

[0078] Figure 1 This is a flow chart of a fast beam shift imaging method provided by one embodiment of the present application. Figure 1 As shown, the fast beam shift imaging method provided by this embodiment may include the following steps:

[0079] S110: Divide the seismic data of the target area into a series of super gathers.

[0080] In the present embodiment, a supergather refers to a combination of shot points and receiver points near a reference point (the supergather center). Because the positions of both shot points and receiver points vary within a supergather, the slope directions of both can be obtained simultaneously. Beam decomposition in a supergather is much easier than in other gathers, so the present invention chooses to perform beam decomposition in a supergather.

[0081] In a specific embodiment, Figure 2 Schematic diagram of 2D and 3D super gather division, such as Figure 2As shown in FIG, dividing the seismic data of the target area into a series of super gathers may include the following steps:

[0082] S111: Shot lines are defined according to the start and end line numbers and the intervals between the vertical and horizontal survey lines of the observation system, and the intersection point on any shot line is used as the center position of the shot point of the super gather. All shot points within the first preset range are selected as the super gather.

[0083] In the embodiment of the present application, the first preset range is a range of 150m in radius centered on the center position of the shot point.

[0084] S112: defining detection point lines according to the vector offset range and the offset sampling interval, taking the intersection of any detection point line as the detection point center position of the super gather, and selecting all detection points within the second preset range as the super gather.

[0085] In the embodiment of the present application, the second preset range is a range of 150 meters around the center of the detection point.

[0086] S120: For each super gather, identifying the arrival time of the local event according to the energy of the local event at the center of the super gather.

[0087] In the embodiments of this application, Figure 3a and Figure 3b As shown, the inventors observed that the coherence surface of cross-events along different slope directions in the super gather is multimodal, and the positions of local extrema accurately correspond to the slope directions of the corresponding events. Therefore, the problem of calculating the slope direction of the cross-event can be transformed into the problem of searching for multiple local extrema of the coherence surface. That is, the coherence of the events along different slope directions can be used as the objective function, which can be transformed into a multimodal optimization problem. This application adopts a variance-based coherence calculation formula as the objective function.

[0088] Figure 4 This is a schematic diagram of identifying the arrival time of local events in a two-dimensional super gather. In a specific embodiment, step S120 may include the following steps:

[0089] S121: Calculate the energy curve of the local event in a single super gather.

[0090] In the embodiment of the present application, the calculation formula of the local event energy curve is:

[0091]

[0092] Where E represents the event energy, t represents the calculation time position, J represents the total number of super gather seismic data, j represents the channel number index, K represents the half-width of the calculation window, k represents the index of the time sampling point in the calculation window, and D represents the seismic data. It is a constant greater than or equal to 1, usually 3, and |*| indicates the absolute value.

[0093] S122: Determine a local maximum value of the energy curve according to the local event energy curve.

[0094] In the embodiment of the present application, the local maximum position of the energy curve can be roughly obtained from Figure 4 It is observed that its exact solution can be calculated by the corresponding software.

[0095] S123: Determine the arrival time of the local event according to the local maximum value of the energy curve.

[0096] In the embodiment of the present application, a local maximum value can be obtained according to the calculation formula of the local event energy curve, and the local maximum value is the arrival time of the local event.

[0097] S130: Determine the slope direction of the local event at the arrival time position of the local event, where the slope direction includes the shot point slope direction and the detection point slope direction.

[0098] In the embodiments of the present application, a three-step multimodal optimization method is used to calculate the slope direction of the local event. The three-step multimodal optimization method specifically includes: (1) a differential evolution algorithm based on neighborhood species; (2) a distance-based individual screening algorithm; and (3) a Hooke-Jeeves pattern search algorithm. This method avoids the over-search problem that occurs when using only the niche method under high-precision computational requirements.

[0099] In a specific embodiment, determining the slope direction of the local event according to the arrival time of the local event may include the following steps:

[0100] S131: The distribution of individuals in the population is driven to the vicinity of the local optimal solution through the differential evolution algorithm based on neighborhood species.

[0101] In an embodiment of the present application, the neighborhood species-based differential evolution algorithm is a multi-peak optimization algorithm that introduces the neighborhood species niche method into the differential evolution algorithm, which can drive population individuals to evolve towards their respective local optimal solutions.

[0102] Specifically, the neighborhood species-based differential evolution algorithm includes:

[0103] Step 1: Randomly initialize and generate N population individuals, and calculate their coherence values.

[0104] In this step, the N population individuals represent N different local event slope directions.

[0105] Step 2: Sort the individuals in the population in ascending order according to the coherence value.

[0106] In this step, the coherence value is calculated using a variance-based coherence calculation formula.

[0107] Step 3: When the sorted population is non-empty,

[0108] The individual with the smallest coherence value among the untreated population individuals is selected as the seed of the species;

[0109] Find the m individuals closest to the seed of the species and set them as one species;

[0110] Delete the processed population individuals in the current population.

[0111] Step 4: Traverse all species, perform mutation and crossover operations, generate test individuals, compare the coherence value of the test individual with the coherence value of the individual closest to it, and if the coherence value of the test individual is smaller, replace the individual closest to it with it.

[0112] Step 5: Merge all species to re-form the population.

[0113] Step 6: If the maximum number of iterations is reached, terminate; otherwise, return to step 2.

[0114] Furthermore, in step 4, the mutation formula is

[0115]

[0116] Among them, V i =[V i,1 ,V i,1 ,…,V i,D ] represents the mutation vector, D represents the dimension of the solution, r1, r2 and r3 are three different random individual indices in the species, and are the same as the current target vector X i The index is different, F represents the scaling factor, which is used to control the amplification effect of the deviation, 0≤F≤2, and is generally taken as 0.9.

[0117] The crossover formula is

[0118]

[0119] Among them, U i,j It represents the component of the jth dimension of the test individual with index i. CR is the control parameter, which is usually 0.1. rand j is a random number uniformly distributed in the range [0,1], and k is a natural number randomly selected between [1,D] to ensure that U i At least from V i Get a portion.

[0120] In step 3, m is generally between 1 / 20 and 1 / 5 of the number of individuals in the population N.

[0121] Repeated experiments found that:

[0122] For 2D supergathers, when the population size is set to 25 and the maximum number of iterations is set to 25 (the number of objective function calls is 625 at this time), the neighborhood species-based differential evolution algorithm can ensure that an excellent solution appears near each local optimal solution.

[0123] For 3D supergathers, when the population size is set to 135 and the maximum number of iterations is set to 60 (the number of objective function calls is 8100), the neighborhood species-based differential evolution algorithm can ensure that an excellent solution appears near each local optimal solution.

[0124] S132: Screening the population individuals near the local optimal solution through the distance-based individual screening algorithm to obtain the excellent individuals closest to each local optimal solution.

[0125] Specifically, the distance-based individual screening algorithm includes the following steps:

[0126] Step 1: Set the maximum threshold C_limit of the coherence of individuals in the population and the minimum threshold D_limit of the distance between individuals in the population.

[0127] Step 2: Traverse all individuals in the population, and for those individuals whose coherence value is greater than C_limit, set their coherence value to 10.0.

[0128] Step 3: Calculate the distance between all population individuals whose coherence value is less than 10.0. If the distance is less than or equal to D_limit, set the coherence value of the individual with the larger coherence value to 10.0.

[0129] Step 4: Output the population individuals whose coherence value is less than 10.0.

[0130] In a specific embodiment, C_limit is generally set to 0.95, and D_limit is generally set to 0.00005.

[0131] S133: The Hooke-Jeeves pattern search algorithm is used to ensure that the selected excellent individuals converge to the local optimal solution with certainty.

[0132] In the embodiments of the present application, the Hooke-Jeeves pattern search algorithm is a local optimization algorithm characterized by not requiring the calculation of the derivative of the objective function during the calculation process. This method ensures that the selected excellent individuals deterministically converge to the local optimal solution, obtaining high-precision calculation results. Each iteration of this method alternates between axial movement and pattern movement. The purpose of axial movement is to find a favorable direction for the objective function to descend, while pattern movement accelerates the search along the favorable direction.

[0133] Figure 5 This is a schematic diagram of the Hooke-Jeeves pattern search. Figure 5 As shown, the Hooke-Jeeves pattern search algorithm provided in this embodiment specifically includes the following steps:

[0134] Step 1: Given an initial point X (1) , n coordinate directions e1,e2,…,e n , initial step size d, acceleration factor α≥1, reduction factor β∈(0,1), error accuracy ε>0, initialize Y (1) =X (1) , k=1, j=1;

[0135] Step 2: If the coherence value C(Y (j) +de j ) <C(Y (j) ), then positive axial movement is performed, that is, Y (j+1) =Y (j) +de j , proceed to step 4; otherwise, proceed to step 3;

[0136] Step 3: If the coherence value C(Y (j) -de j ) <C(Y (j) ), then negative axial movement is performed, that is, Y (j+1) =Y (j) -de j Otherwise, let Y (j+1) =Y (j) ;

[0137] Step 4: If j≤n, set j=j+1 and go to step 2; otherwise, go to step 5;

[0138] Step 5: If C(Y (n+1) ) <C(X (k) ), then proceed to step 6; otherwise, proceed to step 7;

[0139] Step 6: Let X (k+1) =Y (n+1) , perform mode movement, that is, Y (1) =X (k+1)+α(X (k+1) -X (k) ), set k=k+1, j=1, and go to step 2;

[0140] Step 7: If d < ε, stop the iteration and get point X (k) Otherwise, let

[0141] d=βd,Y (1) =X (k) ,X (k+1) =X (k)

[0142] Then let k=k+1, j=1, go to step 2, and repeat the above operation.

[0143] S140: Filter the obtained local optimal solution using a distance-based individual screening algorithm.

[0144] In the embodiment of the present application, the slope direction of the local event finally obtained in step S130 is screened again using an individual screening algorithm based on distance, thereby further improving the calculation accuracy.

[0145] S150: Superimposing the local events into wavelets along the slope direction of the local events to construct a beam.

[0146] In the embodiment of the present application, this step superimposes the local phase axis into a sub-wave along the slope direction of the shot point and the detection point, and simultaneously saves its arrival time, shot point and detection point position (super gather center position), shot point and detection point slope direction to construct the beam.

[0147] S160: Perform kinematic and dynamic ray tracing along the slope direction of the corresponding local event at the center position of the super gather of the beam, and determine the imaging position of the beam according to the time and space positioning conditions.

[0148] In the embodiment of the present application, the purpose of this step is to achieve beam positioning. The specific process of beam positioning is as follows:

[0149] S161: Solve the kinematic ray tracing equation and the dynamic ray tracing equation to obtain the ray paths from the shot point and the receiver point respectively and the complex travel time of the Gaussian beam around the ray path.

[0150] In the embodiment of the present application, the kinematic ray tracing equation is as follows:

[0151]

[0152] Among them, i=1, 2, 3 represents the coordinate direction, x i represents the coordinates of the ray, t represents the travel time along the ray, and p irepresents the slowness component of the ray, and v represents the velocity at the ray position.

[0153] The dynamic ray tracing equation is as follows:

[0154]

[0155] Among them, P, Q and V are all 2×2 symmetric matrices, P and Q are complex-valued dynamic ray parameters, and V represents the second-order derivative of the velocity v with respect to the coordinates of the ray center.

[0156] In an embodiment of the present application, the kinematic ray tracing equation and the dynamic ray tracing equation are solved by the fourth-order Runge-Kutta method, and the complex travel time of the ray path from the shot point, the ray path from the detection point, and the Gaussian beam around the above ray path can be obtained.

[0157] S162: Searching for a point closest to the arrival time of the local event near the nearest position between the ray paths from the shot point and the receiver point, and using it as the imaging position of the beam.

[0158] In the embodiment of the present application, the closest position between two ray paths from the shot point and the receiver point is the spatial positioning condition, which is as follows:

[0159] |S(A)-R(B)| <d max

[0160] Where S(A) is a point on the ray from the shot point, R(B) is a point on the ray from the receiver point, and d max represents the maximum correlation distance.

[0161] Assume that the travel times from the shot point and the receiver point are t s (M) and t r (M), where M is the imaging point position of the beam. The t s (M) and t r (M). The point closest to the arrival time of the local event is the time positioning condition, which is as follows:

[0162] t s (M)+t r (M)≈t event

[0163] Among them, t event represents the arrival time of the local event.

[0164] like Figure 6As shown in FIG, the imaging position of the beam is to find the point M closest to the arrival time of the local phase axis near the closest position (i.e., S(A) and R(B)) of the two rays (i.e., the ray with an asterisk from the shot point and the ray with a triangle from the detection point) when the spatial positioning condition is met.

[0165] S170: At the imaging position of the beam, the beam is expanded according to the real travel time and the virtual travel time of the Gaussian beam at the center position of the super gather to obtain a beam expansion result.

[0166] In the embodiment of the present application, the purpose of this step is to achieve beam spreading.

[0167] First, the calculation formula for Gaussian beam is

[0168]

[0169] Among them, u GB represents the Gaussian beam, x represents the coordinates of the imaging point, x L represents the emission position of the Gaussian beam (shot point or detection point position), ω represents the angular frequency, represents the complex-valued travel time of the Gaussian beam, represents the complex-valued amplitude of the Gaussian beam, and n represents the normal distance from the imaging point to the ray.

[0170] Figure 7 is a schematic diagram of beam expansion, as shown in Figure 7 As shown in Figure 2, the specific process of beam expansion is as follows:

[0171] S171: For any imaging point near the imaging position of the beam, determine the sum of the real part and the sum of the imaginary part of the complex travel time of the Gaussian beam from the shot point and the receiver point.

[0172] In the embodiment of the present application, since P and Q are complex-valued dynamic ray parameters, Has a real part and an imaginary part.

[0173] S172: Using the sum of the real parts of the Gaussian beam, interpolate in the wavelet of the beam to obtain the image field value of the imaging point.

[0174] S173: Attenuate the image field value using the sum of the imaginary parts of the Gaussian beam, thereby obtaining a beam expansion result.

[0175] Through the above steps, beam expansion can be performed on all imaging points near the imaging position of the beam to obtain a beam expansion result.

[0176] S180: Accumulate the beam expansion results of each super gather and output the imaging result.

[0177] In the embodiment of the present application, after all super gathers are processed according to steps S110 to S160, this step may be performed to accumulate the beam expansion results of each super gather and output the imaging result.

[0178] Figure 8 For Marmousi2 speed model.

[0179] Figure 9 This is the migration result diagram of Marmousi2 velocity model. Figure 9 As shown, Figure 9 -a is the Gaussian beam offset, Figure 9 -b is the fast beam offset when all crossed events are considered, Figure 9 -c is the fast beam offset when only a single event is considered, Figure 9 -d is Figure 9 -b with Figure 9 -c difference result.

[0180] Figure 10 This is a partial enlarged view of the Marmousi2 velocity model migration result. Figure 10 middle, Figure 10 -a, Figure 10 -b, Figure 10 -c, Figure 10 -d corresponds to Figure 9 -b; Figure 10 -e, Figure 10 -f, Figure 10 -g, Figure 10 -h corresponds to Figure 9 -c.

[0181] Figure 11 Velocity profiles at different locations of the 3D SEG / EAGE salt model. Figure 11 In a, z = 1020m, Figure 11 -b, x = 4620m, Figure 11 -c, x = 6260m, Figure 11 -d, y = 6280m.

[0182] Figure 12 This is the migration result of the 3D SEG / EAGE salt model at z = 1020m. Figure 12 -a is the Gaussian beam offset, Figure 12 -b is the fast beam offset when all crossed events are considered, Figure 12 -c is the fast beam offset when only a single event is considered, Figure 12 -d is Figure 12 -b with Figure 12 -c difference result.

[0183] Figure 13 This is the migration result of the 3D SEG / EAGE salt model at x = 4620m. Figure 13 -a is the Gaussian beam offset, Figure 13 -b is the fast beam offset when all crossed events are considered, Figure 13 -c only considers the fast beam deviation of a single event, Figure 13 -d is Figure 13 -b with Figure 13 -c difference result.

[0184] Figure 14 This is the migration result of the 3D SEG / EAGE salt model at x = 6260m. Figure 14 -a is the Gaussian beam offset, Figure 14 -b is the fast beam offset when all crossed events are considered, Figure 14 -c is the fast beam offset when only a single event is considered, Figure 14 -d is Figure 14 -b with Figure 14 -c difference result.

[0185] Figure 15 This is the migration result of the 3D SEG / EAGE salt model at y=6280m. Figure 15 -a is the Gaussian beam offset, Figure 15 -b is the fast beam offset when all crossed events are considered, Figure 15 -c is the fast beam offset when only a single event is considered, Figure 15 -d is Figure 15 -b with Figure 15 -c difference result.

[0186] Figure 16 This is a schematic diagram of the structure of a fast beam shift imaging device provided by an embodiment of the present application. The device can be in the form of software and / or hardware. Figure 16 As shown, the fast beam migration imaging device 1600 provided in this embodiment includes: a data processing module 1601, an arrival time identification module 1602, a slope direction determination module 1603, a beam decomposition module 1604, a beam positioning module 1605, a beam expansion module 1606 and an imaging output module 1607, wherein:

[0187] The data processing module 1601 is used to divide the seismic data of the target area into a series of super gathers.

[0188] The arrival time identification module 1602 is used to identify the arrival time of the local event according to the energy of the local event at the center of each super gather.

[0189] The slope direction determining module 1603 is used to determine the slope direction of the local event at the arrival time position of the local event. The slope direction includes the shot point slope direction and the detection point slope direction.

[0190] The beam decomposition module 1604 is configured to superimpose the local events into wavelets along the slope direction of the local events to construct a beam.

[0191] The beam positioning module 1605 is used to perform kinematic and dynamic ray tracing along the slope direction of the corresponding local event at the center position of the super gather of the beam, and determine the imaging position of the beam according to the time and space positioning conditions.

[0192] The beam spreading module 1606 is configured to spread the beam at the imaging position of the beam according to the real travel time and the virtual travel time of the Gaussian beam at the center position of the super gather to obtain a beam spreading result.

[0193] The imaging output module 1607 is used to accumulate the beam expansion results of each super gather and output the imaging results.

[0194] In some embodiments, the slope direction determination module 1603 is specifically configured to:

[0195] The differential evolution algorithm based on neighborhood species is used to drive the distribution of individuals in the population to the vicinity of the local optimal solution;

[0196] The population individuals near the local optimal solution are screened by the distance-based individual screening algorithm to obtain the excellent individuals closest to each local optimal solution;

[0197] The Hooke-Jeeves pattern search algorithm is used to ensure that the selected excellent individuals converge to the local optimal solution with certainty.

[0198] In some embodiments, the fast beam migration imaging device 1600 further includes a screening module for further screening the local optimal solution output by the slope direction determination module 1603 using a distance-based individual screening algorithm.

[0199] In some embodiments, the data processing module 1601 is specifically configured to:

[0200] The shot lines are defined according to the start and end line numbers and the line intervals of the vertical and horizontal survey lines of the observation system. The intersection point on any shot line is taken as the center position of the shot point of the super gather, and all shot points within the first preset range are selected as the super gather.

[0201] The detection point lines are defined according to the vector offset range and the offset sampling interval, and the intersection of any detection point line is used as the detection point center position of the super gather. All detection points within the second preset range are selected as the super gather.

[0202] In some embodiments, the time identification module 1602 is specifically configured to:

[0203] Calculate the energy curve of the local event in a single super gather using the following formula:

[0204]

[0205] Where t represents the calculation time position, J represents the total number of super gather seismic data, j represents the channel number index, K represents the half-width of the calculation window, k represents the index of the time sampling point in the calculation window, and D represents the seismic data. It is a constant greater than or equal to 1. |*| indicates taking the absolute value.

[0206] Determine the local maximum value of the energy curve according to the energy curve of the local event;

[0207] The arrival time of the local event is determined according to the local maximum value of the energy curve.

[0208] In some embodiments, the beam positioning module 1605 is specifically configured to:

[0209] Solve the kinematic ray tracing equation and the dynamic ray tracing equation to obtain the complex travel time of the ray paths from the shot point and the receiver point and the Gaussian beam around the ray paths;

[0210] The point closest to the arrival time of the local event is found near the closest position between the ray paths from the shot point and the receiver point, and is used as the imaging position of the beam.

[0211] In some embodiments, the beam spreading module 1606 is specifically configured to:

[0212] For any imaging point near the imaging position of the beam, determine the sum of the real part and the sum of the imaginary part of the complex travel time of the Gaussian beam from the shot point and the receiver point;

[0213] The image field value of the imaging point is obtained by interpolating the sum of the real parts of the Gaussian beam in the wavelet of the beam;

[0214] The image field value is attenuated using the sum of the imaginary parts of the Gaussian beam to obtain the beam expansion result.

[0215] The fast beam shift imaging device 1600 provided in this embodiment can execute the fast beam shift imaging method provided in the above method embodiment. Its implementation principle and technical effects are similar and will not be described in detail in this embodiment.

[0216] Figure 17 This is a schematic diagram of the structure of a fast beam shift imaging system provided by one embodiment of the present application. Figure 17As shown, the fast beam migration imaging system 1700 provided in this embodiment includes: at least one processor 1701 and a memory 1702. Optionally, the system also includes a communication component 1703. The processor 1701, the memory 1702, and the communication component 1703 are connected via a bus 1704.

[0217] During the specific implementation process, at least one processor 1701 executes the computer-executable instructions stored in the memory 1702, so that the at least one processor 1701 performs the above method.

[0218] The specific implementation process of the processor 1701 can be found in the above method embodiment. Its implementation principle and technical effects are similar and will not be repeated here in this embodiment.

[0219] In the above embodiments, it should be understood that the processor may be a central processing unit (CPU), other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), etc. A general-purpose processor may be a microprocessor or any conventional processor. The steps of the method disclosed in the present invention may be directly implemented by a hardware processor or implemented by a combination of hardware and software modules in the processor.

[0220] The memory may include a high-speed memory (Random Access Memory, RAM), and may also include a non-volatile memory (NVM), such as at least one disk memory.

[0221] The bus can be an Industry Standard Architecture (ISA) bus, a Peripheral Component Interconnect (PCI) bus, or an Extended Industry Standard Architecture (EISA) bus. Buses can be classified into address buses, data buses, and control buses. For ease of illustration, the buses in the drawings of this application are not limited to just one bus or just one type of bus.

[0222] The present application also provides a computer program product, including a computer program, which implements the above method when executed by a processor.

[0223] The present application also provides a computer-readable storage medium, in which computer-executable instructions are stored. When a processor executes the computer-executable instructions, the above method is implemented.

[0224] The above-mentioned readable storage medium can be implemented by any type of volatile or non-volatile memory 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 memory, flash memory, magnetic disk or optical disk. The readable storage medium can be any available medium that can be accessed by a general-purpose or special-purpose computer.

[0225] An exemplary readable storage medium is coupled to a processor so that the processor can read information from the readable storage medium and write information to the readable storage medium. The readable storage medium may also be an integral part of the processor. The processor and the readable storage medium may be located in an application specific integrated circuit (ASIC). Of course, the processor and the readable storage medium may also exist in a device as discrete components.

[0226] The division of units is merely a logical functional division; actual implementations may employ alternative divisions, such as combining or integrating multiple units or components into another system, or omitting or disabling certain features. Furthermore, any direct coupling or communication connection shown or discussed may be an indirect coupling or communication connection between devices or units, either through an interface, electrical, mechanical, or other means.

[0227] Units described as separate components may or may not be physically separate, and components shown as units may or may not be physical units, that is, they may be located in one place or distributed across multiple network units. Some or all of these units may be selected to achieve the purpose of this embodiment according to actual needs.

[0228] In addition, each functional unit in each embodiment of the present invention may be integrated into one processing unit, or each unit may exist physically separately, or two or more units may be integrated into one unit.

[0229] If the function is implemented in the form of a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, or the part that contributes to the prior art, or the part of the technical solution, can be embodied in the form of a software product. The computer software product is stored in a storage medium and includes several instructions for enabling a computer device (which can be a personal computer, server, or network device, etc.) to execute all or part of the steps of the various embodiments of the present invention. The aforementioned storage medium includes: U disk, mobile hard disk, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), disk or optical disk, and other media that can store program code.

[0230] Those skilled in the art will appreciate that all or part of the steps in the above-described method embodiments can be implemented using hardware associated with program instructions. The aforementioned program can be stored in a computer-readable storage medium. When executed, the program performs the steps of the above-described method embodiments. The aforementioned storage medium includes various media capable of storing program code, such as ROM, RAM, magnetic disks, or optical disks.

[0231] Finally, it should be noted that those skilled in the art will readily identify other embodiments of the present invention after considering the specification and practicing the invention disclosed herein. The present invention is intended to cover any variations, uses, or adaptations of the present invention that follow the general principles of the present invention and include common knowledge or customary techniques in the art not disclosed herein. The present invention is not limited to the precise structure described above and illustrated in the accompanying drawings, and various modifications and variations may be made without departing from the scope thereof. The scope of the present invention is limited solely by the appended claims.

Claims

1. A fast beam shift imaging method, characterized in that: include: Divide the seismic data of the target area into a series of super gathers; For each super gather, the arrival time of the local event is identified based on the energy of the local event at the center of the super gather; Determining the slope direction of the local event at the arrival time position of the local event, wherein the slope direction includes the shot point slope direction and the detection point slope direction; The local events are superimposed into wavelets along the slope direction of the local events to construct a beam; Perform kinematic and dynamic ray tracing along the slope direction of the corresponding local event at the center of the beam super gather, and determine the imaging position of the beam according to the temporal and spatial positioning conditions; At the imaging position of the beam, the beam is expanded according to the real travel time and virtual travel time of the Gaussian beam at the center of the super gather to obtain the beam expansion result; Accumulate the beam expansion results of each super gather and output the imaging result; Determining the slope direction of the local event at the arrival time position of the local event includes: The differential evolution algorithm based on neighborhood species is used to drive the distribution of individuals in the population to the vicinity of the local optimal solution; The population individuals near the local optimal solution are screened by the distance-based individual screening algorithm to obtain the excellent individuals closest to each local optimal solution; The Hooke-Jeeves pattern search algorithm is used to ensure that the selected excellent individuals converge to the local optimal solution with certainty.

2. The method according to claim 1, characterized in that Also includes: The obtained local optimal solutions are further screened using a distance-based individual screening algorithm.

3. The method according to claim 1, characterized in that The seismic data of the target area is divided into a series of super gathers, including: Shot lines are defined according to the start and end line numbers and line intervals of the vertical and horizontal survey lines of the observation system. The intersection point on any shot line is used as the center position of the shot point of the super gather. All shot points within the first preset range are selected as the super gather. The detection point lines are defined according to the vector offset range and the offset sampling interval, and the intersection of any detection point line is used as the detection point center position of the super gather. All detection points within the second preset range are selected as the super gather.

4. The method according to claim 1, wherein The method of identifying the arrival time of the local event according to the energy of the local event at the center of each super gather comprises: Calculate the energy curve of the local event in a single super gather; determining a local maximum value of an energy curve according to the energy curve of the local event; The arrival time of the local event is determined according to the local maximum value of the energy curve.

5. The method according to claim 1, wherein The kinematic and dynamic ray tracing is performed along the slope direction of the corresponding local event at the center position of the super gather of the beam, and the imaging position of the beam is determined according to the time and space positioning conditions, including: Solving the kinematic ray tracing equation and the dynamic ray tracing equation to obtain the ray paths from the shot point and the receiver point respectively and the complex travel time of the Gaussian beam around the ray path; The point closest to the arrival time of the local event is found near the closest position between the ray paths from the shot point and the receiver point, and is used as the imaging position of the beam.

6. The method according to claim 5, characterized in that The step of expanding the beam at the imaging position of the beam according to the real travel time and the virtual travel time of the Gaussian beam at the center position of the super gather to obtain a beam expansion result includes: For any imaging point near the imaging position of the beam, determine the sum of the real part and the sum of the imaginary part of the complex travel time of the Gaussian beam from the shot point and the receiver point; The image field value of the imaging point is obtained by interpolating the sum of the real parts of the Gaussian beam in the wavelet of the beam; The image field value is attenuated by utilizing the sum of the imaginary parts of the Gaussian beam, thereby obtaining a beam expansion result.

7. A rapid beam shift imaging device, characterized in that include: A data processing module for dividing the seismic data of the target area into a series of super gathers; An arrival time identification module is used to identify the arrival time of the local event based on the energy of the local event at the center of each super gather; A slope direction determination module is used to determine the slope direction of the local event at the arrival time position of the local event, wherein the slope direction includes the shot point slope direction and the detection point slope direction; A beam decomposition module is used to superimpose local events into wavelets along the slope direction of the local events to construct a beam; A beam positioning module is used to perform kinematic and dynamic ray tracing along the slope direction of the corresponding local event at the center of the beam super gather, and determine the imaging position of the beam according to the temporal and spatial positioning conditions; A beam spreading module is used to spread the beam at the imaging position of the beam according to the real travel time and the virtual travel time of the Gaussian beam at the center position of the super gather to obtain a beam spreading result; An imaging output module is used to accumulate the beam expansion results of each super gather and output the imaging results; The slope direction determination module is specifically configured to: The differential evolution algorithm based on neighborhood species is used to drive the distribution of individuals in the population to the vicinity of the local optimal solution; The population individuals near the local optimal solution are screened by the distance-based individual screening algorithm to obtain the excellent individuals closest to each local optimal solution; The Hooke-Jeeves pattern search algorithm is used to ensure that the selected excellent individuals converge to the local optimal solution with certainty.

8. A rapid beam shift imaging system, characterized in that include: at least one processor, and a memory communicatively coupled to the at least one processor; The memory stores computer-executable instructions; The processor executes the computer-executable instructions stored in the memory to implement the method according to any one of claims 1 to 6.

9. A computer-readable storage medium, characterized in that The computer-readable storage medium stores computer-executable instructions, which are used to implement the method according to any one of claims 1 to 6 when executed by a processor.

Images