Least square reverse time migration method and device for three-dimensional anisotropic medium

Abstract

The invention discloses a least square reverse time migration method and device for a three-dimensional anisotropic medium, and the method comprises the steps: building a forward equation of the three-dimensional anisotropic medium for a seismic block to be subjected to reverse time migration, and obtaining a background wave field under the three-dimensional anisotropic medium through regional forward modeling; determining a three-dimensional anisotropic medium adjoint wave field by using a pre-established wave equation of the adjoint wave field; according to the adjoint wave field and the background wave field, determining an initial gradient by using a pre-established least square reverse time migration cross-correlation equation, and obtaining a migration result according to the initial gradient; and updating an imaging result through a conjugate gradient method according to the offset result to obtain a reflected wave record.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of petroleum geology basic research, and in particular to a least square reverse time migration method and device for three-dimensional anisotropic medium. BACKGROUND

[0002] Reverse time migration is a high-precision migration imaging method, which can adapt to complex structures, can image turning waves, multiple waves, prismatic waves, etc., and has no angle limitation. However, reverse time migration also has its own inherent shortcomings. Since a large number of finite difference calculations are required in reverse time migration, the reverse time migration algorithm has the defect of high intensive calculation.

[0003] Considering that anisotropy exists in underground medium in general, and the current exploration means mainly rely on three-dimensional observation, therefore, the least square reverse time migration under the current common two-dimensional assumption and the three-dimensional isotropic least square reverse time migration are not in line with the actual situation.

[0004] The current anisotropic least square reverse time migration under two-dimensional conditions has certain effect. This method does not consider the influence of the change of reflection interface and anisotropic parameters in the direction perpendicular to the survey line on imaging, and the two-dimensional processing method is extremely low in efficiency and cannot efficiently process three-dimensional data.

[0005] In addition, three-dimensional data processing has large amount of calculation. Considering the application of high-performance calculation based on a graphics processing unit (GPU) to seismic imaging calculation application, the energy consumption and calculation waiting time caused by large amount of calculation are reduced.

[0006] Therefore, a three-dimensional anisotropic VTI medium least square reverse time migration method based on GPU parallel is proposed to solve the current problems. SUMMARY

[0007] The present application provides a least square reverse time migration method and device for three-dimensional anisotropic medium. The method applies the first-order velocity stress equation in three-dimensional anisotropic VTI medium as the basis for wave field continuation, further develops the corresponding three-dimensional GPU accelerated least square reverse time migration algorithm, improves the resolution and signal-to-noise ratio of seismic data migration imaging, and can eliminate lateral reflection and the artifacts caused by the assumption of isotropic medium that cannot simulate the wave field in anisotropic medium.

[0008] In the first aspect, the present application provides a least square reverse time migration method and device for three-dimensional anisotropic medium, which comprises:

[0009] A forward equation of three-dimensional anisotropic medium is established for a seismic block to be subjected to reverse time migration, and a background wave field under three-dimensional anisotropic medium is obtained by regional forward calculation;

[0010] determining the accompanying wave field of the three-dimensional anisotropic medium by using a pre-established wave equation of the accompanying wave field;

[0011] determining an initial gradient by using a pre-established least square reverse time migration cross-correlation equation according to the accompanying wave field and the background wave field, and obtaining a migration result according to the initial gradient;

[0012] updating the imaging result by a conjugate gradient method according to the migration result, and obtaining a reflection record.

[0013] In a second aspect, an embodiment of the present application further provides a least square reverse time migration device for a three-dimensional anisotropic medium, the device comprising: a memory and a processor; the memory is configured to store a program for performing the least square reverse time migration for the three-dimensional anisotropic medium, and the processor is configured to read and execute the program for performing the least square reverse time migration for the three-dimensional anisotropic medium, and execute the method in any one of the above embodiments.

[0014] In a third aspect, an embodiment of the present application further provides a computer readable storage medium, and the computer readable storage medium stores a data processing program, and the data processing program is executed by a processor to perform the least square reverse time migration method for the three-dimensional anisotropic medium in any one of the above embodiments.

[0015] Compared with the related art, the present application provides a least square reverse time migration method and device for a three-dimensional anisotropic medium, the method comprising: establishing a forward equation for a three-dimensional anisotropic medium for a seismic block to be subjected to reverse time migration, obtaining a background wave field under the three-dimensional anisotropic medium by regional forward calculation; determining the accompanying wave field of the three-dimensional anisotropic medium by using a pre-established wave equation of the accompanying wave field; determining an initial gradient by using a pre-established least square reverse time migration cross-correlation equation according to the accompanying wave field and the background wave field, and obtaining a migration result according to the initial gradient; and updating the imaging result by a conjugate gradient method according to the migration result, and obtaining a reflection record. The first-order velocity stress equation in the three-dimensional anisotropic medium is applied as the basis for wave field continuation, and a corresponding three-dimensional GPU accelerated least square reverse time migration algorithm is further developed, so that the resolution and signal-to-noise ratio of the seismic data migration imaging are improved, and the artifacts caused by the assumption of lateral reflection and isotropic assumption that cannot simulate the wave field in the anisotropic medium are eliminated.

[0016] Other features and advantages of the present application will be set forth in the following description, and in part will become apparent to those skilled in the art from the following description, or can be learned by practice of the present application. Other advantages of the present application can be realized and obtained by means of the schemes described in the specification and the drawings. BRIEF DESCRIPTION OF DRAWINGS

[0017] The accompanying drawings are used to provide an understanding of the technical solutions of the present application, and constitute a part of the specification, and are used together with the embodiments of the present application to explain the technical solutions of the present application, and do not constitute a limitation on the technical solutions of the present application.

[0018] Figure 1 A flow chart of a least square reverse time migration method for a three-dimensional anisotropic medium of embodiments of the present application;

[0019] Figure 2 A schematic diagram of a least square reverse time migration device for a three-dimensional anisotropic medium of embodiments of the present application;

[0020] Figure 3 A schematic diagram of a three-dimensional velocity model in some example embodiments;

[0021] Figure 4 A three-dimensional anisotropic parameter epsilon model in some example embodiments;

[0022] Figure 5 A three-dimensional anisotropic parameter delta model in some example embodiments;

[0023] Figure 6 A three-dimensional isotropic reverse time migration result in some example embodiments;

[0024] Figure 7 A three-dimensional isotropic least square reverse time migration result in some example embodiments;

[0025] Figure 8 A three-dimensional VTI reverse time migration result in some example embodiments;

[0026] Figure 9 A three-dimensional VTI least square reverse time migration result in some example embodiments;

[0027] Figure 10 A partial display diagram of a three-dimensional VTI reverse time migration result in some example embodiments;

[0028] Figure 11 A partial display diagram of a three-dimensional VTI least square reverse time migration result in some example embodiments;

[0029] Figure 12 A three-dimensional velocity field in some example embodiments;

[0030] Figure 13 A three-dimensional isotropic reverse time migration result in some example embodiments;

[0031] Figure 14 A three-dimensional isotropic least square reverse time migration result in some example embodiments;

[0032] Figure 15 Three-dimensional VTI reverse-time migration results in some example embodiments;

[0033] Figure 16 Three-dimensional VTI least-squares reverse-time migration results in some example embodiments;

[0034] Figure 17 Partial display of three-dimensional VTI reverse-time migration results in some example embodiments;

[0035] Figure 18 Partial display of three-dimensional VTI least-squares reverse-time migration results in some example embodiments;

[0036] Figure 19 GPU parallel-based three-dimensional VTI medium least-squares reverse-time migration method flowchart in some example embodiments. DETAILED DESCRIPTION

[0037] The present application describes a number of embodiments, but the description is exemplary rather than limiting and it will be apparent to those of ordinary skill in the art that numerous more embodiments and implementations are possible within the scope of the embodiments described in the present application. Although a number of possible combinations of features have been set forth herein, many more combinations are possible. Unless specifically intended otherwise, any feature or element of any embodiment can be used in combination with any other feature or element of any other embodiment, or in combination with any other feature or element of the same embodiment.

[0038] The present application includes and contemplates combinations of features and elements known to those of ordinary skill in the art. The embodiments, features and elements disclosed herein can also be combined with any conventional features or elements to form unique inventive solutions that are within the scope of the claims. Any feature or element of any embodiment can also be combined with features or elements from other inventive solutions to form another unique inventive solution that is within the scope of the claims. Therefore, it should be understood that any feature shown and / or discussed in the present application can be used, alone or in any appropriate combination. Accordingly, the embodiments are not to be restricted, except as by the appended claims and their equivalents. Furthermore, various modifications and changes can be made within the scope of the appended claims.

[0039] Furthermore, in describing representative embodiments, the specification can have presented the method and / or process as a particular sequence of steps. However, to the extent that the method or process depends on the performance of certain steps, the steps need not be performed in the order described. Other sequences of steps can be possible, and are contemplated by those of ordinary skill in the art. Therefore, the particular order in which the steps are presented is not limiting. Furthermore, for the sake of clarity, the description has not attempted to include all the possible combinations of implementations, procedures, or steps that can be implemented or performed by one of ordinary skill in the art. It is understood that the scope of the description is intended to encompass all such possible combinations.

[0040] At present, the least square reverse time migration technology mainly develops two-dimensional isotropic least square reverse time migration, two-dimensional anisotropic least square reverse time migration and three-dimensional isotropic least square reverse time migration. The least square reverse time migration is mainly for optimizing the imaging effect of the reverse time migration method, and this method has certain effect in actual application. However, the basic least square reverse time migration imaging method is usually established on the assumption of isotropic homogeneous acoustic medium, and such assumption can basically meet the requirement of basic imaging and has been widely applied in the industry. However, the actual underground medium is complex and changeable. The actual underground medium usually has multiple complex properties such as elasticity, viscosity, anisotropy and inhomogeneity. Therefore, the conventional assumption does not meet the actual situation of the underground, and when these properties have a strong influence on the seismic wave dynamics and kinematics characteristics, the traditional simple acoustic assumption will cause a series of negative effects such as imaging artifact and noise. Therefore, many experts and scholars studying LSRTM also try to further develop the least square reverse time migration imaging LSRTM under complex medium based on the seismic wave propagation mechanism of complex medium.

[0041] The embodiment of the present application provides a least square reverse time migration method of three-dimensional anisotropic medium, as shown in the figure, the method comprises steps S100-S130: Figure 1

[0042] Step S100: establishing a forward equation of three-dimensional anisotropic medium for a seismic block to be subjected to reverse time migration, and obtaining a background wave field under three-dimensional anisotropic medium through regional forward;

[0043] Step S110: determining a three-dimensional anisotropic medium accompanying wave field by using a wave equation of the accompanying wave field established in advance;

[0044] Step S120: determining an initial gradient by using a least square reverse time migration cross-correlation equation established in advance according to the accompanying wave field and the background wave field, and obtaining a migration result according to the initial gradient; ​

[0045] Step S130: updating the imaging result by the conjugate gradient method according to the offset result to obtain the current iteration of the reflected wave record.

[0046] In the embodiment, the GPU processor under the control of the CPU thread is used to perform forward modeling of single-shot seismic data corresponding thereto in the three-dimensional VTI medium to obtain the background wave field in the three-dimensional VTI medium.

[0047] In the embodiment, one CPU thread controls one GPU processor; and one GPU processor processes one single-shot seismic data.

[0048] In an example embodiment, the initial gradient is determined according to the accompanying wave field and the background wave field by using a pre-established least square reverse time migration cross-correlation equation, and the offset result is obtained according to the initial gradient, which includes:

[0049] The single-shot cross-correlation result under each GPU is determined according to the accompanying wave field and the background wave field by using a least square reverse time migration cross-correlation equation.

[0050] The initial gradient of the current iteration is determined according to the single-shot cross-correlation results under all GPUs.

[0051] The operation results in the multiple GPUs are merged by an MPI reduction function according to the initial gradient to obtain the offset result.

[0052] In the embodiment, the accompanying wave field in the three-dimensional anisotropic medium is determined by using a pre-established wave equation of the accompanying wave field to perform reverse time continuation of the shot record residual; and the GPU and the single-shot record are one-to-one corresponding, that is, one CPU thread controls one GPU, and one GPU processes one single-shot data.

[0053] In the embodiment, the LSRTM cross-correlation equation in the three-dimensional acoustic VTI medium is derived in combination with the accompanying state method, then the reverse time continuation of the shot record residual is performed by using the wave equation of the accompanying wave field, the obtained accompanying wave field and the background field are cross-correlated, and finally the single-shot cross-correlation results calculated in each GPU are reduced and integrated to determine the initial gradient of the current iteration.

[0054] In an example embodiment, the wave equation of the background wave field in the three-dimensional anisotropic medium is:

[0055]

[0056] In the equation, u0, v0 and w0 are background intermediate variables, p0 and q0 are scalar background wave fields, s acts as a source to promote the forward propagation of the wave field, t is time, x, y and z represent three coordinate axes in a three-dimensional orthogonal coordinate system, V0 is background velocity, p is density, and ε and δ are Thomson parameters.

[0057] In an example embodiment, the wave equation for the wavefield accompanying the three-dimensional anisotropic medium is:

[0058]

[0059] In the above equation, u adj , v adj , w adj represent the velocity components in three different directions in the accompanying wavefield; p adj , q adj represent two components of the scalar wavefield accompanying the wavefield; p s , q s represent two simulated reflection wavefield records; p obs , q obs represent two actually observed shot records.

[0060] In an example embodiment, the initial gradient is:

[0061]

[0062] In the above equation, V0 is the background velocity, J is the objective function; m is the reflection coefficient, p0 and q0 are two scalar background wavefield components;

[0063] In the above equation, J is the objective function:

[0064]

[0065] J is the objective function, p obs , q obs represent two actually observed shot records, C represents the calculation of new reflection wave records, d represents the actually observed reflection wave record data.

[0066] In an example embodiment, the imaging result is updated by the conjugate gradient method according to the migration result to obtain the reflection wave record of the current iteration, including:

[0067] The imaging result is updated by the conjugate gradient method according to the migration result; the updated imaging result is taken as the reflection coefficient, and the equation of the reverse migration processing is used for reflection wave simulation;

[0068] The reflection wave record of the current iteration is determined according to the reflection wave simulation result.

[0069] In an example embodiment, the equation of the reverse migration processing is:

[0070]

[0071] In the above equation, u s, v s , w s The disturbance quantity of the velocity three-component is represented by Δv, and the disturbance quantity of the velocity field is represented by ΔV.

[0072] In an example embodiment, determining the reflection wave record of the current iteration according to the updated imaging result comprises:

[0073] Performing reverse migration processing on the updated imaging result as the reflection coefficient;

[0074] Obtaining the reflection wave record of the current iteration according to the reverse migration processing.

[0075] In an example embodiment, the GPU parallel-based three-dimensional VTI medium least square reverse time migration method process comprises:

[0076] Step S10: performing forward modeling of the corresponding region of the three-dimensional VTI medium by using the GPU under the control of the CPU to obtain the background wave field under the three-dimensional VTI medium;

[0077] Step S11: establishing the wave equation of the accompanying wave field under the three-dimensional VTI medium of the corresponding region;

[0078] Step S12: establishing the objective function and determining the least square reverse time migration cross-correlation equation under the three-dimensional VTI medium;

[0079] Step S13: performing reverse time continuation of the shot record residual by using the wave equation of the accompanying wave field to determine the accompanying wave field;

[0080] Step S14: determining the single-shot cross-correlation result under each GPU by using the least square reverse time migration cross-correlation equation according to the accompanying wave field and the background wave field;

[0081] Step S15: determining the initial gradient of the current iteration according to the single-shot cross-correlation results under all GPUs;

[0082] Step S16: obtaining the updated imaging result according to the initial gradient and the imaging result of the last iteration;

[0083] Step S17: determining the reflection wave record of the current iteration according to the updated imaging result;

[0084] Step S18: determining whether the iteration end condition is met according to the residual of the reflection wave record of the current iteration and the actual observed record, if not, re-executing steps 13-17; if yes, outputting the migration result.

[0085] In a second aspect, the embodiments of the present application further provide a device for least square reverse time migration of three-dimensional anisotropic medium, characterized in that the device comprises a memory and a processor; the memory is configured to store a program for least square reverse time migration of three-dimensional anisotropic medium, and the processor is configured to read and execute the program for least square reverse time migration of three-dimensional anisotropic medium, and execute the method of any one of the above embodiments.

[0086] In a third aspect, the embodiments of the present application further provide a computer readable storage medium, which stores a data processing program, and the data processing program is executed by a processor to execute the least square reverse time migration method of three-dimensional anisotropic medium of any one of the above embodiments.

[0087] Example one

[0088] As Figure 19 shown, the flow of the GPU parallel-based three-dimensional VTI medium least square reverse time migration method is specifically described as follows:

[0089] Step 1. Parameter field input;

[0090] The related parameters such as migration velocity and anisotropic parameter field are input.

[0091] Step 2. Thread division according to the total number of GPU;

[0092] Single-shot migration task is assigned to one thread in CPU, one CPU thread controls one GPU, and one GPU processes one single-shot migration data.

[0093] Step 3. Cutting processing of the calculation region according to the observation system;

[0094] According to the range covered by the detector in each single shot, the three-dimensional data body is cut, the calculation range of single shot is compressed, and the cut data body is transmitted into each GPU.

[0095] Step 4. Parameter field extension;

[0096] When wave field simulation is performed, the calculated region has a boundary, which will cause false reflection at the boundary that does not conform to the actual situation. Therefore, before calculation, a number of absorbing layers are added to the periphery of the cut parameter field to absorb wave field energy and prevent false reflection. Because of the added absorbing layer, the parameter field in the layer is empty, so the value of the outermost parameter field is assigned to the absorbing layer to make the absorbing layer also have a parameter field.

[0097] Step 5. Input local parameter field;

[0098] Step 6. Initialization of the array for storing wave field;

[0099] The wavefield is calculated, and a space is opened in the computer to store this wavefield, which is an array.

[0100] Initialization of the array storing the wavefield means that the values in the array are all set to 0.

[0101] Step 7. Forward modeling under a three-dimensional VTI medium;

[0102] Step 71. Forward propagation of the wavefield in the background velocity field;

[0103] Step 72. Store the background wavefield at the boundary at each time;

[0104] Step 8. Adjoint source back propagation and cross-correlation imaging under a three-dimensional VTI medium; the adjoint source back propagation under a three-dimensional VTI medium is to use the wave equation of the adjoint wavefield to perform reverse time extension of the shot record residual;

[0105] Step 81. Background wavefield reconstruction using the stored boundary wavefield;

[0106] Step 82. Residual as an adjoint source back propagation, calculate the adjoint wavefield;

[0107] Step 83. Background wavefield and adjoint wavefield cross-correlation to obtain single-shot initial gradient

[0108]

[0109] Step 9. Initial gradient calculation;

[0110] Integrate the single-shot migration results on each GPU using the MPI reduction function to form the total initial gradient.

[0111] Step 10. Update using the conjugate gradient method combined with the initial gradient to obtain the imaging result of this iteration;

[0112] Step 11. Local migration results are transmitted back to the corresponding position in space;

[0113] Step 12. According to the updated imaging result, calculate the new reflected wavefield using Born forward modeling (inverse migration), and calculate the residual between the new reflected wave record and the original observation record;

[0114] Step 13. The residual is passed into step 8 as an adjoint source, and steps 8-12 are repeated until the residual is less than a threshold value or a certain number of iterations is reached, and then the iteration is stopped.

[0115] Example two

[0116] Based on the reverse-time migration principle in the field of seismic exploration, combined with the least square theory, the least square reverse-time migration technology is formed. Considering the limitations of conventional isotropic method and two-dimensional anisotropic method in processing effect, the first-order velocity stress equation in VTI medium is applied as the basis of wave field continuation in the application, and the corresponding three-dimensional GPU accelerated least square reverse-time migration algorithm is further developed, so that the resolution and signal-to-noise ratio of imaging results are improved, and the artifacts caused by lateral reflection and isotropic assumption that cannot simulate wave field in anisotropic medium are eliminated. Based on the three-dimensional VTI medium least square reverse-time migration technology of GPU parallel, the implementation process is as follows:

[0117] Step 1, according to the range covered by the detector in each shot, the three-dimensional data body is cut, and the cut data body is respectively transmitted into each GPU.

[0118] The cut multi-shot task is assigned to multiple CPU threads for simultaneous calculation, wherein each CPU thread controls a GPU, and a GPU processes a single-shot seismic data.

[0119] The conventional CPU parallel code is calculated based on CPU and the running memory of the computer, the CPU single-core computing capacity is strong, but the number of CPU cores is limited, which limits the computing efficiency. The CUDA is parallel calculated based on GPU and the GPU memory, the computing capacity of the GPU core is relatively simple, but the number of GPU cores is much larger than that of the conventional CPU. Unlike the calculation mode of each computing core in CPU parallel bearing all the calculation tasks of grid points in a shot, the calculation mode of GPU is that each computing core bears the calculation task of a grid point in the model. Due to the large number of GPU cores, the computing efficiency is obviously improved compared with CPU parallel.

[0120] In the application, not only the GPU is used for single-shot processing, but also the multi-GPU cooperative parallel algorithm is developed. The basic framework of this algorithm is that the multi-shot task is assigned to multiple CPU threads for simultaneous calculation, wherein each CPU thread controls a GPU, and the GPU bears the main calculation task. After the calculation is completed, the results calculated by multiple GPUs are combined by using the reduction function of MPI, and the calculation results are updated.

[0121] Step 2, the forward simulation of the corresponding region of three-dimensional VTI medium is carried out in each GPU to obtain the boundary wave field information at each time sampling point.

[0122] If the exploration area has obvious anisotropy, RTM or LSRTM must use the differential operator containing anisotropy to accurately simulate the actual propagation of the wave field. Since the forward background wave field is needed in subsequent imaging, the background wave field in the three-dimensional VTI medium needs to be simulated in this step, and its boundary values are saved for subsequent reconstruction.

[0123] Under the pseudo-acoustic wave assumption, the first-order velocity-stress equation in the VTI medium, i.e., formula 1, represents the background acoustic wave equation in the VTI medium, which can be expressed as:

[0124]

[0125] where u0, v0 and w0 are background intermediate variables, p0 and q0 are scalar background wave fields, s acts as a source to promote the forward propagation of the wave field, t is time, and x, y and z represent three coordinate axes in a three-dimensional rectangular coordinate system. The above equation mainly contains four constants: V0 is the background velocity, the density p, the Thomson parameter e and d.

[0126] By discretizing the above equation by finite difference, the background wave field in the three-dimensional VTI medium can be calculated.

[0127] Step 3, calculate the three-dimensional VTI accompanying wave field at the corresponding position in each GPU by inverse time, and reconstruct the forward background wave field, and perform cross-correlation imaging of the two to obtain the initial gradient. After obtaining the initial gradient, the corresponding cross-correlation results in all CPU threads are combined into one by using the MPI reduction function to form the total cross-correlation result.

[0128] The wave equation of the accompanying wave field under three-dimensional VTI is:

[0129]

[0130] Formula 15 obtains the accompanying wave field. According to the results of formula 15 and the background wave field obtained by formula 1, the initial gradient is obtained by using formula 16 to cross-correlate. Through the obtained initial gradient, the corresponding cross-correlation results in all CPU threads are combined into one by using the MPI reduction function to form the total cross-correlation result.

[0131] Step 4, according to the basic theory of the conjugate gradient method, the cross-correlation result obtained in the above step is used to update the imaging result.

[0132] According to the cross-correlation result obtained in step 3, i.e., the initial gradient, the initial gradient obtained is transformed to obtain the update direction; then multiplied by the iteration step, and then added to the imaging result obtained by the last iteration to obtain the imaging result of this time.

[0133] Step 5, after the updated imaging results are taken as the reflection coefficients and subjected to the cutting process similar to that in Step 1, the corresponding GPU is returned to the Born forward corresponding to the position, that is, the de-migration process, to obtain new reflection records, and the specific calculation is as follows:

[0134] According to the Born approximation theory, the wave field, the velocity field and the anisotropy parameter field are considered to be linearly decomposable. The wave field is divided into background wave field and perturbation wave field; the velocity field is divided into background velocity field and perturbation velocity field; the parameters of anisotropic medium are divided into background parameters and perturbation parameters.

[0135]

[0136] wherein p, q are the total fields of two scalar wave fields; V represents the true velocity; ε, δ are the background fields of the true anisotropy parameters.

[0137] p0 and q0 are scalar background wave fields, V0 is the background velocity, ΔV is the perturbation of V, Δε is the perturbation of ε, and Δδ is the perturbation of δ.

[0138] The influence of anisotropy perturbation on the wave field is small, so the perturbation can be ignored, and only the perturbations of the wave field and the velocity are retained. In order to simplify the derivation, the mathematical approximation

[0139]

[0140] Then formula 1 can be equivalent to

[0141]

[0142] wherein u s , v s , w s represent the perturbations of the velocity three components. Similarly, p s and q s represent the perturbations of the wave field.

[0143] Subtracting formula 4 from formula 1, the de-migration equation in three-dimensional VTI medium is obtained as follows:

[0144]

[0145] The new reflection record can be obtained by using the de-migration equation of formula 5.

[0146] Step 6, the residual error between the reflection record and the actual observation record is obtained;

[0147] According to the objective function of LSRTM, the residual error between the reflection record and the actual observation record is obtained:

[0148]

[0149] where C represents the new reflectivity record calculated in step 5, and d represents the actual observed reflectivity record data.

[0150] Step 7, taking the residual as the source input to step 3 as the source of the adjoint equation, repeating the process of step 3-step 6 until the number of iterations is met or the convergence condition (residual less than or equal to the set value) is reached, obtaining the final three-dimensional VTI least squares reverse-time migration processing result.

[0151] In this example, the determination process of the initial gradient is as follows:

[0152] This method is to realize LSRTM in the data domain, so the RMT process of the data residual is approximately equal to the initial gradient of each iteration.

[0153] The migration operator is the inverse operator of the forward operator in mathematics, which is difficult to realize in calculation. According to the adjoint state method, the adjoint of the forward operator can be regarded as the migration operator. Based on the adjoint state method, it can be obtained

[0154] <A adj U adj ,U>=<U adj ,AU> (Equation 6)

[0155] where A adj is the adjoint matrix of the positive operator A. U adj represents the adjoint wave field, and A can be expanded as

[0156]

[0157] Substituting equation 7 into equation 6, the adjoint matrix of the positive operator can be derived:

[0158]

[0159] The difference between LSRTM and RTM is that LSRTM introduces a target function to measure the imaging quality. The purpose of this target function is to find the scattering potential field to minimize the residual between the Born approximation de-migration obtained reflectivity record and the observed data. This makes LSRTM an optimization problem, and the target function is:

[0160]

[0161] where C represents the reflectivity record obtained by de-migration, and d represents the actual observed data. In addition, the scattering potential is defined as

[0162]

[0163] To solve the optimal solution of m, the derivative of the objective function with respect to the gradient can be expressed as

[0164]

[0165] s s is the scattering source.

[0166] Similarly, according to the adjoint state method, the gradient expression can be expressed as

[0167]

[0168] The adjoint wave field can be regarded as the product of the inverse of the adjoint matrix and the source. For the adjoint wave field, the adjoint source is the residual of the observed data, and the adjoint source can be expressed as

[0169] u adj = (A adj ) -1 (C-d) (Formula 13)

[0170] Formula 13 can be further expressed as

[0171] A adj u adj = C-d (Formula 14)

[0172] The adjoint forward operator is brought into Formula 14

[0173]

[0174] Based on the least squares reverse-time migration theory, the gradient expression is:

[0175]

[0176] Using the method of the present example, based on the designed velocity model, and using the same source wavelet to excite at different positions on the ground surface corresponding to the velocity model, the shot record is forward modeled, and the observed record is used for inversion.

[0177] Figure 3 The three-dimensional velocity model of the example has 301 sampling points in the X direction, 301 sampling points in the Y direction, and 301 sampling points in the Z direction, and the sampling interval in the three directions is 11m. Figure 4 The three-dimensional anisotropic parameter epsilon model of the example has 301 sampling points in the X direction, 301 sampling points in the Y direction, and 301 sampling points in the Z direction, and the sampling interval in the three directions is 11m. Figure 5 The three-dimensional anisotropic parameter delta model of the example has 301 sampling points in the X direction, 301 sampling points in the Y direction, and 301 sampling points in the Z direction, and the sampling interval in the three directions is 11m. Figures 3 to 5The true velocity field and anisotropy field used for forward modeling.

[0178] The forward modeling record uses a 25Hz Ricker wavelet, an observation time of 2.5s, and a time sampling interval of 1ms. Figure 6 is a three-dimensional isotropic reverse-time migration result, and it can be seen from the migration result that the imaging result cannot correctly reflect the position of the actual reflection surface underground, and the imaging result is basically incorrect. Figure 7 is a three-dimensional isotropic least squares reverse-time migration result, and due to a large number of false images introduced by the isotropic assumption, the least squares reverse-time migration result cannot converge at all, and the imaging result is basically the same as the reverse-time migration result. Figure 8 is a three-dimensional reverse-time migration result based on the VTI assumption, and by comparing the position of the reflection surface in Figure 3 , it can be seen that Figure 8 the result basically reflects the correct reflection surface position, but there are still problems such as insufficient focusing of reflection energy and uneven longitudinal energy. Figure 9 The result obtained by the three-dimensional VTI medium least squares reverse-time migration used in the application is basically similar to Figure 8 , but it is more balanced in energy at the middle and deep reflection positions, and better inverts the underground reflection structure, proving the effectiveness of the application. Figure 10 and Figure 11 are respectively Figure 8 and Figure 9 are partial enlarged displays of and

[0179] , and the partial enlarged view further embodies the high resolution and high energy focusing provided by the method of the application.

[0179] In order to further verify the effectiveness and processing effect of the method of the application, a model combination in which the velocity model and the anisotropy model are inconsistent in structure is used for testing, and in the example, the observation system is moving, and a moving cutting parameter field is also used in the test, and a strategy of distributing to multiple GPUs for calculation is adopted. Figure 12 is the velocity model in the example, the global anisotropy parameter epsilon is 0.25, and the global anisotropy parameter delta is 0.2. The test excites 120 shots, and the observation system is an observation system with 10201 receivers uniformly distributed around the shot point as the center. The total observation time is 3s, and the time sampling interval is 1ms. Figure 13 is the three-dimensional isotropic reverse-time migration result obtained under the model, and it can be seen that there are still a large number of false images around the true reflection surface position, especially in the shallow layer, and this phenomenon is particularly obvious. Figure 14 is a three-dimensional isotropic least squares reverse-time migration result, and similarly, under the condition that the medium assumption is incorrect, the least squares reverse-time migration technology cannot bring obvious improvement to the imaging result. Figure 15It is a three-dimensional VTI reverse-time migration result, and it can be seen that, in combination with the wave field propagation theory in the VTI medium, the imaging algorithm can better correct the false image around the reflection surface. Figure 16 It is the result obtained by the present application, and it can be seen that the present application method not only greatly eliminates the false image in imaging, but also realizes and to some extent improves the amplitude balance and the resolution of the imaging result. Figure 17 and Figure 18 It is a slice comparison of three-dimensional VTI reverse-time migration and three-dimensional VTI least square reverse-time migration, and by comparing the part pointed by the arrow in the two figures, the effectiveness of the present application method is further verified.

[0180] The three-dimensional VTI medium least square reverse-time migration technology based on GPU parallelism has the following technical effects:

[0181] I. The least square theory is used to optimize the reverse-time migration method, and on this basis, the wave field propagation theory in the VTI medium is used to correct the deviation of the least square reverse-time migration under the isotropic assumption. In addition, since the method is developed based on three-dimensional observation, the result is more reasonable compared with the two-dimensional method, and the false image caused by lateral reflection can be effectively eliminated.

[0182] II. Wide range of adaptation: the method is developed for three-dimensional observation system, and can be directly used for processing of three-dimensional data.

[0183] III. Simple operation and easy implementation. Through algorithm integration and external integration of input parameters, only the inversion parameters need to be adjusted in a single interface, and the source code does not need to be opened for modification.

[0184] Those of ordinary skill in the art will realize and understand that all or some of the steps in the methods disclosed above and the functional modules / units in the systems and devices can be implemented as software, firmware, hardware, and appropriate combinations thereof. In hardware implementation, the division between the functional modules / units mentioned in the above description does not necessarily correspond to the division of physical components; for example, one physical component can have multiple functions, or one function or step can be performed by several physical components in cooperation. Some or all of the components can be implemented as software executed by a processor, such as a digital signal processor or a microprocessor, or as hardware, or as an integrated circuit, such as an application-specific integrated circuit. Such software can be distributed on computer-readable media, which can include computer storage media (or non-transitory media) and communication media (or transitory media). As is well known to those of ordinary skill in the art, the term computer storage media includes volatile and non-volatile, removable and non-removable media implemented in any method or technology for storage of information such as computer readable instructions, data structures, program modules or other data. Computer storage media includes, but is not limited to, RAM, ROM, EEPROM, flash memory or other memory technology, CD-ROM, digital versatile disks (DVD) or other optical disk storage, magnetic cassettes, magnetic tapes, magnetic disk storage or other magnetic storage devices, or any other medium which can be used to store the desired information and which can be accessed by a computer. Furthermore, it is common and well understood by those of ordinary skill in the art that communication media typically embodies computer readable instructions, data structures, program modules or other data in a modulated data signal such as a carrier wave or other transport mechanism and can include any information delivery media.

Claims

1. A least-squares reverse-time migration method for three-dimensional anisotropic media, characterized in that, The method includes: A forward modeling equation for a three-dimensional anisotropic medium is established for the seismic block to be reversed in time migration, and the background wave field under the three-dimensional anisotropic medium is obtained by regional forward modeling. The associated wave field in a three-dimensional anisotropic medium is determined by using a pre-established wave equation for the associated wave field. Based on the adjoint wave field and the background wave field, the initial gradient is determined using the pre-established least squares inverse time migration cross-correlation equation, and the migration result is obtained based on the initial gradient. The imaging results are updated using the conjugate gradient method based on the offset results to obtain the reflected wave record.

2. The least-squares reverse-time migration method for three-dimensional anisotropic media according to claim 1, characterized in that, The process involves establishing forward modeling equations for a three-dimensional anisotropic medium in the seismic block to be reverse-time migrated, and obtaining the background wave field under the three-dimensional anisotropic medium through regional forward modeling, including: The GPU processor, controlled by CPU threads, is used to perform forward modeling of the corresponding single-shot seismic data in a three-dimensional VTI medium to obtain the background wave field in a three-dimensional anisotropic medium. One CPU thread controls one GPU processor; one GPU processor handles one seismic data point.

3. The least-squares reverse-time migration method for three-dimensional anisotropic media according to claim 2, wherein determining the initial gradient based on the adjoint wave field and the background wave field using a pre-established least-squares reverse-time migration cross-correlation equation, and obtaining the migration result based on the initial gradient, includes: Based on the accompanying wave field and the background wave field, the single-shot cross-correlation result under each GPU is determined using the least squares inverse time-shift cross-correlation equation. The initial gradient for this iteration is determined based on the inter-camera results across all GPUs. Based on the initial gradient, the computation results from multiple GPUs are merged using the MPI reduction function to obtain the offset result; Among them, the adjoint wave field of the three-dimensional anisotropic medium is determined by the inverse time extension of the shot record residual using the pre-established wave equation of the adjoint wave field.

4. The least-squares reverse-time migration method for three-dimensional anisotropic media according to claim 1, wherein the wave equation of the background wave field in the three-dimensional anisotropic media is: in, u0, v0, and w0 are background intermediate variables, p0 and q0 are scalar background wave fields, s acts as a source that promotes the forward propagation of the wave field, t is time, x, y, and z represent the three coordinate axes in a three-dimensional Cartesian coordinate system, V0 is the background velocity, ρ is the density, and ε and δ are Thomson parameters.

5. The least-squares reverse-time migration method for three-dimensional anisotropic media according to claim 4, characterized in that, The wave equation for the adjoint wave field in the three-dimensional anisotropic medium is as follows: In the above formula, u adj v adj w adj p represents the velocity components in three different directions in the accompanying wave field. adj q adj p represents the two components of the scalar adjoint wave field; s q s This represents the simulated record of the reflected wave field; p obs q obs This refers to the actual gun recordings observed.

6. The least-squares reverse-time migration method for three-dimensional anisotropic media according to claim 5, characterized in that, The initial gradient is: In the above formula, V0 is the background velocity; J is the objective function; m is the reflection coefficient; and p0 and q0 are two scalar background wave field components.

7. The least-squares reverse-time migration method for three-dimensional anisotropic media according to claim 6, characterized in that, The step of updating the imaging results using the conjugate gradient method based on the offset results to obtain the reflected wave record includes: The imaging results are updated using the conjugate gradient method based on the offset results. The updated imaging results are used as reflection coefficients, and the reflected waves are simulated using the equations of the inverse migration process. The reflected wave record is determined based on the simulation results of the reflected wave.

8. The least-squares reverse-time migration method for three-dimensional anisotropic media according to claim 7, wherein the equation for the inverse migration process is: In the above formula, u s v s w s ΔV represents the perturbation of the velocity component, and ΔV represents the perturbation of the velocity field.

9. A least-squares reverse-time migration device for a three-dimensional anisotropic medium, characterized in that, The apparatus includes a memory and a processor; the memory is used to store a program for performing least-squares reverse-time migration of a three-dimensional anisotropic medium, and the processor is used to read and execute the program for performing least-squares reverse-time migration of a three-dimensional anisotropic medium, and to perform the method according to any one of claims 1-8.

10. A computer-readable storage medium storing a data processing program, the data processing program being executed by a processor using the least-squares reverse-time migration method for a three-dimensional anisotropic medium according to any one of claims 1-8.

Images