Fast beam migration imaging method, device and equipment

Abstract

The invention provides a fast beam migration imaging method, device and equipment. The method comprises the following steps: dividing seismic data of a target area into a series of super gathers; for each super gather, identifying the arrival time of the local event according to the energy of the local event; determining the slope direction of the local event at the arrival time position of the local event; overlapping the local events into wavelets along the slope direction of the local events, and constructing a wave beam; kinematics and dynamics ray tracing is carried out at the super gather center position of the wave beam along the slope direction of the corresponding local event, and the imaging position of the wave beam is determined; at the imaging position of the wave beam, the wave beam is expanded according to the real travel time and the virtual travel time of the Gaussian wave beam from the central position of the super gather; and accumulating the beam expansion result of each super gather and outputting an imaging result. According to the method, the imaging effect of a complex structure area can be improved under the condition that the calculation problem of the cross event slope is solved.

Description

Technical Field

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

[0002] The fast beam migration imaging method is to decompose the seismic record into a series of beams (wavelets with specific arrival time, shot and receiver position, and shot and receiver slope direction), and then only need to perform imaging along the determined beam direction (shot and receiver slope direction), and the beam decomposition step only needs to be performed once (independent of speed), which has the characteristics of fast speed and high precision. Fast beam migration mainly includes three steps: beam decomposition, beam positioning and beam expansion. The accurate decomposition of the beam is a prerequisite for the accurate positioning and expansion of the beam, and plays a decisive role in the final imaging quality. In the beam decomposition step, the most critical part is the determination of the local event slope direction (shot and receiver slope direction). The accuracy of the local event slope direction will directly affect 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 multi-dimensional local tilt stacking method to obtain it directly in the super gather; the third is to use the plane wave decomposition filter method to obtain it in the super gather.

[0004] None of the above technologies can handle the situation of event crossing well. However, for seismic data, the phenomenon of event crossing is very common, especially for seismic data from complex structural areas (such as unconformities or faults). Due to 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 calculation efficiency of beam decomposition. Summary of the invention

[0005] In view of the problem that the prior art cannot handle the situation of crossed phase axes well and the beam decomposition calculation efficiency is low, resulting in poor imaging effect 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 according to 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 sub-waves along the slope direction of the local events to construct a beam;

[0011] At the center of the super gather of the beam, kinematic and dynamic ray tracing is performed along the slope direction of the corresponding local event axis, and the imaging position of the beam is determined according to the time and space positioning conditions;

[0012] 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 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 population individuals to the vicinity of the local optimal solution;

[0017] The population individuals near the local optimal solution are screened through 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] The shot lines are defined according to the start and end line numbers and the 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 the shot points within the first preset range are selected as the super gather;

[0023] The detection point line is defined according to the vector offset range and the offset sampling interval, the intersection of any detection point line is taken as the detection point center position of the super gather, and all the 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 the super gather for 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 axis at the center position of the super gather of the beam to determine the imaging position of the beam according to the time and space positioning conditions includes:

[0030] The kinematic ray tracing equation and the dynamic ray tracing equation are solved to obtain the complex travel time of the ray paths from the shot point and the detector point and the Gaussian beam around the ray paths;

[0031] The point closest to the arrival time of the local event axis is found near the nearest 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 using the sum of the imaginary parts of the Gaussian beam, thereby obtaining a beam expansion result.

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

[0037] A data processing module, used 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 axis according to the energy of the local event axis at the center of each super gather;

[0039] A slope direction determination module, 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, used for superimposing the local event axes into sub-waves along the slope direction of the local event axes to construct a beam;

[0041] A beam positioning module is used to perform kinematic and dynamic ray tracing along the slope direction of the corresponding local event axis 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;

[0042] A beam expansion module is used to expand 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;

[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 also 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, because the coherence surface of the cross-phase axis 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 phase axis, therefore, the present application can transform the calculation problem of the slope of the cross-phase axis into the problem of searching for multiple local extreme values ​​of the coherence surface. That is, the coherence of the phase axis along different slope directions can be used as the objective function and transformed into a multi-peak optimization problem. The present application solves the slope of the local phase axis in the fast beam migration beam decomposition link through the differential evolution algorithm based on neighborhood species, the individual screening algorithm based on distance and the Hooke-Jeeves pattern search algorithm, so as to obtain high-precision calculation results. This method can reduce the impact of high-dimensional calculation problems on calculation efficiency in the case of three-dimensional super gathers, so as to improve the imaging effect of complex structural areas when dealing with the calculation problem of the slope of the cross-phase axis. 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 flow chart of a fast beam shift imaging method provided by one embodiment of the present application;

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

[0052] Figure 3a It is a schematic diagram of two-dimensional super gather cross-event seismic data;

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

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

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

[0056] Figure 6 This is a schematic diagram of beam positioning;

[0057] Figure 7 Schematic diagram of beam expansion;

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

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

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

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

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

[0063] Fig.13 This is the migration result diagram of the 3D SEG / EAGE salt model at x=4620m;

[0064] Fig.14 This is the migration result diagram of the 3D SEG / EAGE salt model at x=6260m;

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

[0066] Fig.16 A schematic diagram of the structure of a fast beam shift imaging device provided by an embodiment of the present application;

[0067] Fig.17 A schematic diagram of the structure of a fast beam shift imaging system provided in one embodiment of the present application.

[0068] The above drawings have shown clear embodiments of the present application, which will be described in more detail later. These drawings and text descriptions are not intended to limit the scope of the present application in any way, but to illustrate the concept of the present application to those skilled in the art by referring to specific embodiments. DETAILED DESCRIPTION

[0069] Exemplary embodiments will be described in detail herein, examples of which are shown in the accompanying drawings. When the following description refers to the drawings, the same numbers in different drawings represent the same or similar elements unless otherwise indicated. The implementations described in the following exemplary embodiments do not represent all implementations consistent with the present application. Instead, they are merely examples of devices and methods consistent with some 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 the same or similar items with substantially the same functions and effects. Those skilled in the art can understand that words such as "first" and "second" do not limit the quantity and execution order, and words such as "first" and "second" do not necessarily limit the difference.

[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 as "exemplary" or "for example" in the present application should not be interpreted as being more preferred or more advantageous than other embodiments or designs. 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 "more than one" refers to two or more.

[0072] It should be noted that the “at…” in the embodiments of the present application may be the instant when a certain situation occurs, or may be a period of time after the occurrence of a certain situation, and the embodiments of the present application do not specifically limit this.

[0073] Fast beam migration decomposes seismic records into a series of beams (wavelets with specific arrival times, shot and receiver positions, and shot and receiver slope directions) through beam decomposition, and then only needs to be imaged along the determined beam direction (shot and receiver slope direction). The beam decomposition step only needs to be performed once (independent of speed), which is fast and accurate. Fast beam migration can be two orders of magnitude faster in computational efficiency than the Kirchhoff depth migration commonly used in industrial production.

[0074] Fast beam migration mainly includes 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. In the beam decomposition step, the most critical thing is to obtain the local event axis slope direction (the slope direction of the shot point and the detector point). The accuracy of the local event axis slope direction will directly affect the accuracy of beam decomposition.

[0075] In the prior art, there are three main ways to obtain the local event axis slope direction in the fast beam offset 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 multi-dimensional local tilt stacking method to obtain it directly in the super gather; the third is to use the plane wave decomposition filter method to obtain it in the super gather.

[0076] However, among the above technologies, the first method is affected by the large sampling interval and event axis crossing phenomenon along the transverse line direction 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 ) for uniform sampling will cause serious high-dimensional computing problems. For example, if the number of sampling points along each direction is set to 100, then 10 8 The third method also needs to face the problem of high-dimensional calculation in three-dimensional case, and the plane wave decomposition filter method is difficult to deal with 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 calculation 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. The method can effectively reduce the impact of high-dimensional calculation problems on calculation efficiency in the case of three-dimensional super gathers through the fast beam migration imaging method, and improve the imaging effect of complex structural areas when dealing with the calculation problem of cross-phase axis slope.

[0078] Figure 1 A flow chart of a fast beam shift imaging method provided by an 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 embodiment of the present application, a super gather refers to a combination of shot points and detector points near a reference point (the center position of the super gather). Since the positions of the shot points and the detector points are variable in the super gather, the slope directions of the shot points and the detector points can be obtained at the same time. It is much easier to perform beam decomposition in a super gather than in other gathers, so the present invention chooses to perform beam decomposition in a super gather.

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

[0082] S111: defining shot lines according to the start and end line numbers and the intervals between the vertical and horizontal survey lines of the observation system, taking the intersection point on any shot line as the center position of the shot point of the super gather, and selecting all shot points within the first preset range 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 centered at 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 the present application, Figure 3a and Figure 3b As shown, the inventors observed that the coherence surface of the cross-events along different slope directions in the super gather is multi-peaked, and the position of the local extrema accurately corresponds to the slope direction of the corresponding event. Therefore, the calculation problem of the slope direction of the cross-event can be transformed into a problem of searching for multiple local extrema of the coherence surface, that is, the coherence of the event along different slope directions can be used as the objective function and transformed into a multi-peak optimization problem. The present application adopts the coherence calculation formula based on variance 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 time window, k represents the index of the time sampling point in the calculation time 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 corresponding software.

[0095] S123: Determine the arrival time of the local event axis 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 a calculation formula of a local event energy curve, and the local maximum value is the arrival time of the local event.

[0097] S130: determining 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 embodiment of the present application, a three-step multi-peak optimization method is used to calculate the slope direction of the local phase axis. The three-step multi-peak optimization method specifically includes: (1) a differential evolution algorithm based on neighborhood species; (2) an individual screening algorithm based on distance; (3) a Hooke-Jeeves pattern search algorithm. The above method avoids the problem of excessive search when only the niche method is used under high-precision computing 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 a 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 differential evolution algorithm based on neighborhood species 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.

[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 indexes in the species, and they 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 0.9.

[0117] The crossover formula is

[0118]

[0119] Among them, U i,j represents the component of the jth dimension of the test individual with index i, CR is the control parameter, which is generally 0.1, and 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 serving.

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

[0121] Repeated experiments found that:

[0122] For 2D super gathers, 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 super gathers, 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 at this time), 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 values ​​are 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 values ​​are 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 embodiment of the present application, the Hooke-Jeeves pattern search algorithm is a local optimization algorithm, which is characterized by not needing to calculate the derivative of the objective function during the calculation process. This method can ensure that the selected excellent individuals converge to the local optimal solution with certainty, and obtain 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, and pattern movement is to accelerate the search along the favorable direction.

[0133] Figure 5 Figure 1 is a schematic diagram of the Hooke-Jeeves mode 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) , to move the mode, 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: Screening 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 axis finally obtained in step S130 is screened again by using an individual screening algorithm based on distance, so as to further improve the calculation accuracy.

[0145] S150: Superimpose 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 event axes into sub-waves along the slope directions of the shot and detection points, and simultaneously saves their arrival times, shot and detection point positions (center positions of super gathers), and shot and detection point slope directions to construct a beam.

[0147] S160: Perform kinematic and dynamic ray tracing along the slope direction of the corresponding local event axis 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, and 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 detector 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 the 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-valued 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: Find a point that is closest to the arrival time of the local event axis near the nearest position between the ray paths from the shot point and the detection point, and use 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 detection 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. Through the paraxial rays, t s (M) and t r (M). The point closest to the arrival time of the local phase axis 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. 1 , 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 rays with asterisks from the shot point and the rays with triangles from the detection point) when the spatial positioning conditions are 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 the Gaussian beam is

[0168]

[0169] Among them, u GB represents the Gaussian beam, x represents the coordinates of the imaging point, and 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, such as 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 detection 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 to output an 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] Fig. 9 This is the migration result diagram of the Marmousi2 velocity model. Fig. 9 As shown, Fig. 9 -a is the Gaussian beam offset, Fig. 9 -b is the fast beam deviation when all crossed events are considered, Fig. 9 -c is the fast beam deviation when only a single event axis is considered, Fig. 9 -d is Fig. 9 -b with Fig. 9 -c difference result.

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

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

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

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

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

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

[0186] Fig.16 This is a schematic diagram of the structure of a fast beam shift imaging device provided by an embodiment of the present application, and the device may be in the form of software and / or hardware. Fig.16 As shown, the fast beam shift 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, and the slope direction includes the shot point slope direction and the detection point slope direction.

[0190] The beam decomposition module 1604 is used to superimpose the local events into sub-waves 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 axis 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 expansion module 1606 is used to expand 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.

[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 population individuals to the vicinity of the local optimal solution;

[0196] The population individuals near the local optimal solution are screened through 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 shift 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 used to:

[0200] The shot line is 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 the shot points within the first preset range are selected as the super gather.

[0201] The detection point line is defined according to the vector offset range and the offset sampling interval, the intersection of any detection point line is taken as the detection point center position of the super gather, and all the 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 axis in a single super gather, and the calculation formula is:

[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 calculation time window half width, k represents the index of the time sampling point in the calculation time window, and D represents the seismic data. is a constant greater than or equal to 1, and |*| indicates taking the absolute value.

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

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

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

[0209] The kinematic ray tracing equation and the dynamic ray tracing equation are solved 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 axis is found near the nearest 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, and its implementation principle and technical effect are similar, which will not be described in detail in this embodiment.

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

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

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

[0219] In the above embodiments, it should be understood that the processor may be a central processing unit (CPU), or 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 invention may be directly implemented as being executed by a hardware processor, or may be executed 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 (Non-volatile Memory, NVM), such as at least one disk memory.

[0221] The bus may be an Industry Standard Architecture (ISA) bus, a Peripheral Component Interconnect (PCI) bus, or an Extended Industry Standard Architecture (EISA) bus, etc. The bus may be divided into an address bus, a data bus, a control bus, etc. For ease of representation, the bus in the drawings of the present application is not limited to only one bus or 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 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 memory, flash memory, magnetic disk or optical disk. The readable storage medium can be any available medium that can be accessed by a general 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 can also be a component of the processor. The processor and the readable storage medium can be located in an application specific integrated circuit (Application Specific Integrated Circuits, referred to as: ASIC). Of course, the processor and the readable storage medium can also exist in the device as discrete components.

[0226] The division of units is only a logical function division, and there may be other divisions in actual implementation, such as multiple units or components can be combined or integrated into another system, or some features can be ignored or not executed. Another point is that the mutual coupling or direct coupling or communication connection shown or discussed can be an indirect coupling or communication connection through some interface, device or unit, which can be electrical, mechanical or other forms.

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

[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, including several instructions for a computer device (which can be a personal computer, server, or network device, etc.) to perform all or part of the steps of the methods of each embodiment 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, etc. Various media that can store program codes.

[0230] Those skilled in the art can understand that all or part of the steps of implementing the above-mentioned method embodiments can be completed by hardware related to program instructions. The aforementioned program can be stored in a computer-readable storage medium. When the program is executed, the steps of the above-mentioned method embodiments are executed; and the aforementioned storage medium includes: ROM, RAM, disk or optical disk and other media that can store program codes.

[0231] Finally, it should be noted that those skilled in the art will readily conceive of 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, which follow the general principles of the present invention and include common knowledge or customary technical means in the art not disclosed by the present invention, are not limited to the precise structure described above and shown in the drawings, and may be modified and changed in various ways without departing from the scope thereof. The scope of the present invention is limited only 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 according to 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 sub-waves along the slope direction of the local events to construct a beam; At the center of the super gather of the beam, kinematic and dynamic ray tracing is performed along the slope direction of the corresponding local event axis, and the imaging position of the beam is determined according to the time and space positioning conditions; 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 the beam expansion result; The beam expansion results of each super gather are accumulated to output the imaging result.

2. The method according to claim 1, characterized in that 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 population individuals to the vicinity of the local optimal solution; The population individuals near the local optimal solution are screened through 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.

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

4. 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: The shot lines are defined according to the start and end line numbers and the 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 the shot points within the first preset range are selected as the super gather; The detection point line is defined according to the vector offset range and the offset sampling interval, the intersection of any detection point line is taken as the detection point center position of the super gather, and all the detection points within the second preset range are selected as the super gather.

5. The method according to claim 1, characterized in that The method of identifying the arrival time of the local event according to the energy of the local event at the center of the super gather for 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.

6. The method according to claim 1, characterized in that The kinematic and dynamic ray tracing is performed along the slope direction of the corresponding local event axis 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: The kinematic ray tracing equation and the dynamic ray tracing equation are solved to obtain the complex travel time of the ray paths from the shot point and the detector point and the Gaussian beam around the ray paths; The point closest to the arrival time of the local event axis is found near the nearest position between the ray paths from the shot point and the receiver point, and is used as the imaging position of the beam.

7. The method according to claim 6, 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 the 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 using the sum of the imaginary parts of the Gaussian beam, thereby obtaining a beam expansion result.

8. A rapid beam shift imaging device, characterized in that: include: A data processing module, used 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 axis according to the energy of the local event axis at the center of each super gather; A slope direction determination module, 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, used for superimposing the local event axes into sub-waves along the slope direction of the local event axes 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 axis 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; A beam expansion module is used to expand 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; The imaging output module is used to accumulate the beam expansion results of each super gather and output the imaging results.

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

10. 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 7 when executed by a processor.

Images