Multiple reverse time migration method based on convolution interference principle

Abstract

The invention provides a multiple reverse time migration method based on a convolution type interference principle, and the method comprises the steps: 1, inputting a speed field, and building an observation system; step 2, generating a primary wave which does not contain a direct wave by adopting a surface absorption boundary condition; step 3, performing primary seismic interference convolution processing on the primary waves to obtain first-order multiples; and step 4, performing conventional multiple imaging and hierarchical multiple imaging by using the obtained first-order multiple. The multiple reverse time migration method based on the convolution interference principle is different from the previous multiple separation technology, and the multiple is predicted by directly using the primary wave. In order to prove the correctness of the method, the predicted multiples are compared with shot records obtained by adopting free earth surface boundary conditions. The multiple reverse time migration is compared with the conventional reverse time migration.

Description

Technical Field

[0001] The invention relates to the technical field of petroleum exploration, and in particular to a multiple wave reverse time migration method based on the convolution interference principle. Background Art

[0002] Conventional reverse time migration is mainly for primary waves with strong energy, rather than multiple waves. Multiple waves are usually removed as a kind of noise. However, multiple waves have a long propagation path and contain richer geological information. Multiple wave reverse time migration has good imaging effect on special geological conditions such as subsalt structures. However, the previous multiple wave construction technology mainly separated multiple waves from the shot records containing primary waves and multiple waves. Therefore, there is currently a lack of technology to directly predict multiple waves with primary waves.

[0003] In the Chinese patent application with application number: CN201510846550.X, a method and system for processing seismic data by prestack multiple wave reverse time migration is involved, and the method includes the following steps: Step 1: Acquire seismic data; Step 2: Start multiple hadoop data nodes on the cluster, divide the seismic data into multiple data blocks, and send each data block to a data node on the cluster; Step 3: Perform prestack multiple wave reverse time migration calculation on the seismic data on each data node to obtain a temporary imaging result, and store the temporary imaging result in a distributed file system; Step 4: Merge the temporary imaging results in the distributed file system to obtain a common imaging point gather as an imaging result; Step 5: Output the imaging result to generate an imaging profile. The invention realizes the prestack multiple wave reverse time migration calculation of seismic data by adopting the Mapreduce framework and GPU acceleration. The method is not only convenient for processing large-scale and ultra-large-scale seismic data, but also can improve the efficiency of migration processing. However, the invention uses prestack multiple wave migration imaging instead of converting multiple waves into primary waves, which is obviously different from the present invention and cannot solve the technical problem that the present invention wants to solve.

[0004] In the Chinese patent application with application number: CN201410188065.3, a reverse time migration method for three-dimensional VSP data is involved, which is a method for reverse time migration imaging of three-dimensional VSP data based on the three-dimensional acoustic wave equation. The present invention can make full use of the rich wave field information in the VSP data and the imaging ability of the reverse time migration method for complex wave fields, realize accurate imaging of complex structures around the well, and provide favorable guidance for the exploration of complex oil and gas reservoirs. The method fully considers the different wave field information actually existing underground (multiple waves, diffraction waves, etc.) and does not simplify the wave equation. The method uses optimized finite difference coefficients to solve the partial differential equation by central difference, which has higher accuracy and more stable calculation than the traditional method. The application of this method to the processing of three-dimensional VSP seismic data is of great significance to the imaging of small faults near the well, reservoir prediction, micro-structure implementation, and the exploration and development of concealed oil and gas reservoirs. However, the invention performs reverse time migration imaging on the VSP observation system and does not take targeted measures for imaging multiple waves, which may result in illusions. This is obviously different from the present invention and cannot solve the technical problem that the present invention is intended to solve.

[0005] In the Chinese patent application with application number: CN202210709847.1, a seismic wavefield migration imaging method of interlayer multiple Fourier finite difference is involved. It is improved based on the full wavefield simulation technology of wavefield extension. At the layer position of each stratum, the wavefields of interlayer multiple waves of different orders are simulated and obtained respectively. According to the propagation path of interlayer multiple waves of different orders, the interlayer multiple waves of different orders are back-transmitted, and the wavefields of interlayer multiple waves of different orders are cross-correlated imaged in the upgoing wavefield and the downgoing wavefield respectively. The input data includes geological model data and seismic records containing only interlayer multiple waves. Through wavefield extension, setting of earthquake source, artificial boundary reflection suppression, and cyclic recursive wavefield order simulation, the seismic wavefield migration imaging of interlayer multiple waves of Fourier finite difference is realized. The invention converts the interlayer multiple wave seismic signal, which is used as a noise signal in the existing imaging method, into an effective signal, performs seismic wavefield migration imaging, and can be used to restore underground structural information. This invention utilizes wave field separation technology for imaging and is based on the Fourier finite difference migration method, which is obviously different from the present invention and cannot solve the technical problem that the present invention is intended to solve.

[0006] In the Chinese patent application with application number: CN201610798688.1, a GPU reverse time migration method and system based on multi-stream optimization is involved, which relates to the field of geophysical exploration. The present invention relates to the following steps: acquiring seismic data; performing reverse time imaging processing on the seismic data, and calculating the data loading kernel function, wavelet back propagation kernel function and data back propagation kernel function in the processing process by a first-level concurrent calculation method; calculating the imaging condition kernel function and the noise filter kernel function in the reverse time imaging process of the seismic data by a second-level concurrent calculation method; obtaining the imaging result according to the calculation result, and generating the imaging profile. The pre-stack reverse time migration calculation of seismic data is accelerated by adopting a multi-stream parallel optimization strategy. This method can not only process ultra-large-scale seismic data, but also improve the efficiency of migration processing. The invention uses the back propagation kernel function for reverse time migration imaging, but is not aimed at multiple wave imaging, which is obviously different from the present invention and cannot solve the technical problem that the present invention wants to solve.

[0007] The above existing technologies are all significantly different from the present invention and fail to solve the technical problem we want to solve. Therefore, we have invented a new multiple wave reverse time migration method based on the convolution interference principle. Summary of the invention

[0008] The purpose of the present invention is to provide a multiple wave reverse time migration method based on the convolution interference principle which has a reasonable design, overcomes the shortcomings of the prior art and has good effects.

[0009] The object of the present invention can be achieved by the following technical measures: a multiple wave reverse time migration method based on the convolution interference principle, the multiple wave reverse time migration method based on the convolution interference principle comprises:

[0010] Step 1, input velocity field and establish observation system;

[0011] Step 2, using the surface absorbing boundary condition to generate a primary wave without direct waves;

[0012] Step 3, performing a seismic interferometric convolution process on the primary wave to obtain a first-order multiple wave;

[0013] Step 4, using the obtained first-order multiple waves to perform conventional multiple wave imaging and fractional multiple wave imaging.

[0014] The purpose of the present invention can also be achieved by the following technical measures:

[0015] In step 1, the source and receiver points are placed on the surface for excitation, and the shot spacing and trace spacing are set to three grid points.

[0016] In step 3, in 3D space, for the acoustic wave medium, the acoustic wave Green's function in the frequency domain satisfies the three-dimensional Helmholtz equation:

[0017]

[0018] in, represents the gradient operator, k represents the wave number, ω represents the angular frequency, δ represents the step function, s represents the source, g represents the detection point, and G represents the Green's function;

[0019] Formula (1) is a second-order partial differential equation with two solutions, namely the causal Green's function G(g|s,ω) and the non-causal Green's function G(g|s,ω) * ; The positions of the excitation point and the detection point in the Green's function can be interchanged, that is, G(g|s,ω)=G(s|g,ω), which conforms to the reciprocity theorem.

[0020] In step 3, G(x|A,ω) represents the Green's function of the seismic wave response excited by point A inside the medium volume V and received by point x inside the volume V. According to formula (1), the Helmholtz equation of points A and B inside the medium can be defined, namely:

[0021]

[0022]

[0023] in, represents the gradient operator, k represents the wave number, ω represents the angular frequency, δ represents the step function, G and G0 represent Green's functions, x, A and B represent the coordinate positions; wave number k = ω / v(g), k0 = ω / v0(g), x∈V, since the medium is isotropic and has uniform velocity, if A and B are inside the medium, then v(x) = v0(x), and if they are outside, then v(x)≠v0(x).

[0024] In step 3, multiply both sides of formula (2) by G0(x|B,ω), and multiply both sides of formula (3) by G(x|A,ω) and then subtract them to obtain:

[0025]

[0026] in, represents the gradient operator, ω represents the angular frequency, δ represents the step function, G and G0 represent the Green's function, x, A and B represent the coordinate position; since the point x is inside the volume V, the method of integration by parts can be obtained:

[0027]

[0028]

[0029] in, represents the gradient operator, ω represents the angular frequency, G and G0 represent the Green's function, and x, A and B represent the coordinate positions.

[0030] In step 3, substitute formula (5) and formula (6) into formula (4) to obtain:

[0031]

[0032] in, represents the gradient operator, ω represents the angular frequency, G and G0 represent the Green's function, x, A and B represent the coordinate position;

[0033] Using Gauss's theorem to integrate the closed surface S of volume V, we can get:

[0034]

[0035] in, represents the gradient operator, ω represents the angular frequency, G and G0 represent the Green's function, x, A and B represent the coordinate position, Represents the unit external normal vector of the boundary of surface S.

[0036] In step 3, the far-field approximation theory is introduced. The incident wave in the homogeneous medium will cause the scatterer to scatter, which can be expressed in the far field as:

[0037] G(x|B,ω)~m(φ,θ)e ikr / r (9)

[0038] Where ω represents the angular frequency, G represents the Green's function, x and B represent the coordinate positions, r represents the distance from B to the scattering point x in the scatterer, and m(φ,θ) is only a function of the angular coordinates and therefore has no effect on the partial derivative of r. The far field approximation of the gradient term in formula (10) is:

[0039]

[0040] in, represents the gradient operator, ω represents the angular frequency, G represents the Green's function, x and A represent the coordinate position, represents the unit external normal vector of the boundary of the surface S, r represents the distance from the scatterer A to the scattering point x, and k represents the wave number.

[0041] In step 3, because there are:

[0042]

[0043] Where ω represents the angular frequency, G represents the Green's function, x and A represent the coordinate positions, r represents the distance from A inside the scatterer to the scattering point x far away, and k represents the wave number;

[0044] So we have:

[0045]

[0046] Where ω represents the angular frequency, G represents the Green's function, x and A represent the coordinate positions, r represents the distance from A inside the scatterer to the scattering point x far away, and k represents the wave number. represents the unit external normal vector of the surface S boundary, represents the gradient operator, i represents the imaginary unit;

[0047] Similarly, there are:

[0048]

[0049] Where ω represents the angular frequency, G0 represents the Green's function, x and B represent the coordinate positions, r represents the distance from B inside the scatterer to the scattering point x far away, and k represents the wave number. represents the unit external normal vector of the surface S boundary, represents the gradient operator and i represents the imaginary unit.

[0050] In step 3, substitute formula (12) (13) into formula (8) to obtain:

[0051]

[0052] Among them, i represents the imaginary unit, G represents the Green's function, k represents the wave number, ω represents the angular frequency, x, A and B represent the coordinate points, S represents the integral curve, and * represents taking the conjugate.

[0053] In step 4, the steps of multiple-order RTM are:

[0054] Step 41. Since the original shot record of the free surface boundary obtains a mixed wave of primary waves and multiple waves, the original shot record of the surface absorption boundary is used to obtain only a primary wave for subsequent processing operations;

[0055] Step 42. Based on the primary wave record, apply convolution type seismic interferometry to predict the first order multiple waves, predict the second order multiple waves based on the predicted first order multiple waves, and deduce the N order multiple waves from the first order;

[0056] Step 43. Based on the zero-order multiple waves and the first-order multiple waves, migration calculation is performed as the forward wave field and the reverse wave field to obtain the first-order multiple wave reverse time migration imaging section;

[0057] Step 44. Use k-1 order multiple waves and k order multiple waves as forward wave field and reverse wave field to perform migration calculation respectively until k=N, and then obtain reverse time migration imaging sections of multiple waves of each order.

[0058] In step 4, the formula for reverse time migration of multiple waves of different orders and the imaging conditions of multiple waves of different orders can be expressed as:

[0059]

[0060] Wherein, m represents the order of multiple waves, F represents the forward wave field, B represents the reverse wave field, M represents the wave field value, x and z represent the coordinate points, and t represents the time. Here, only the first-order multiple waves are used, so N=1.

[0061] The multiple-wave reverse time migration method based on the convolution-type interference principle in the present invention comprises the following steps: inputting a velocity field and establishing an observation system; using the surface absorption boundary condition to generate a primary wave without a direct wave; performing a seismic interference convolution process on the primary wave to obtain a first-order multiple wave; and using the obtained first-order multiple wave for conventional multiple wave imaging and order-multiple wave imaging. This method is different from the previous multiple wave separation technology, but directly predicts multiple waves with a single wave. In order to prove the correctness of the method, the predicted multiple waves are compared with the shot records obtained using the free surface boundary condition. The multiple wave reverse time migration is compared with the conventional reverse time migration. The multiple wave reverse time migration method based on the convolution-type interference principle in the present invention uses a single wave to realize the prediction of multiple waves, and the multiple wave reverse time migration improves the imaging accuracy of the subsalt structure. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 It is a schematic diagram of a three-layer horizontal model in a specific embodiment of the present invention;

[0063] Figure 2 It is a schematic diagram of a single shot record of a three-layer horizontal model in a specific embodiment of the present invention;

[0064] Figure 3 It is a schematic diagram of reverse time migration of a three-layer horizontal model in a specific embodiment of the present invention;

[0065] Figure 4 A schematic diagram of an inclined layer model in a specific embodiment of the present invention;

[0066] Figure 5 It is a schematic diagram of a single shot record of a diagonal layer model in a specific embodiment of the present invention;

[0067] Figure 6 It is a schematic diagram of reverse time migration of a slant layer model in a specific embodiment of the present invention;

[0068] Figure 7 A schematic diagram of a salt dome model in a specific embodiment of the present invention;

[0069] Figure 8 It is a schematic diagram of a single shot record of a salt dome model in a specific embodiment of the present invention;

[0070] Fig. 9 It is a schematic diagram of reverse time migration of a salt dome model in a specific embodiment of the present invention;

[0071] Fig.10 The flowchart is a specific embodiment of the multiple wave reverse time migration method based on the convolution interference principle of the present invention. DETAILED DESCRIPTION

[0072] It should be noted that the following detailed descriptions are exemplary and are intended to provide further explanation of the present invention. Unless otherwise specified, all technical and scientific terms used herein have the same meanings as those commonly understood by those skilled in the art to which the present invention belongs.

[0073] It should be noted that the terms used herein are only for describing specific embodiments and are not intended to limit exemplary embodiments according to the present invention. As used herein, unless the context clearly indicates otherwise, the singular form is also intended to include the plural form. In addition, it should be understood that when the terms "comprising" and / or "including" are used in this specification, it indicates the presence of features, steps, operations and / or combinations thereof.

[0074] like Fig.10 As shown, Fig.10 The flowchart of the multiple wave reverse time migration method based on the convolution interference principle of the present invention is as follows. The multiple wave reverse time migration method based on the convolution interference principle includes:

[0075] Step 101: Input the velocity field and establish an observation system;

[0076] Step 102: Generate a primary wave without a direct wave using a surface absorbing boundary condition;

[0077] Step 103: Perform a seismic interferometric convolution process on the primary wave to obtain a first-order multiple wave;

[0078] Step 104: Use the obtained first-order multiple waves to perform conventional multiple wave imaging and fractional-order multiple wave imaging.

[0079] The following are several specific embodiments of the present invention.

[0080] Example 1

[0081] In a specific embodiment 1 of the present invention, the multiple wave reverse time migration method based on the convolution interference principle includes the following steps:

[0082] Step 1: Input velocity field and establish observation system;

[0083] The source and receiver points are placed on the surface for excitation, and the shot spacing and trace spacing are set to three grid points.

[0084] Step 2: Use the surface absorbing boundary condition to generate a primary wave without direct waves;

[0085] Step 3: Perform a seismic interferometric convolution process on the primary wave to obtain the first-order multiple waves;

[0086] In 3D space, for acoustic wave media, the acoustic wave Green's function in the frequency domain satisfies the three-dimensional Helmholtz equation:

[0087]

[0088] in, represents the gradient operator, k represents the wave number, ω represents the angular frequency, δ represents the step function, s represents the source, g represents the detection point, and G represents the Green's function.

[0089] Formula (1) is a second-order partial differential equation with two solutions, namely the causal Green's function G(g|s,ω) and the non-causal Green's function G(g|s,ω) * . The positions of the excitation point and the detection point in the Green function can be interchanged, that is, G(g|s,ω)=G(s|g,ω), which conforms to the reciprocity theorem. The reciprocity theorem is a theory that states that the wave field has a reciprocal property, that is, the excitation point and the receiving point of the acoustic wave field are interchanged, and the wave field remains unchanged. The acoustic wave reciprocity theorem includes two parts: the convolution type reciprocity theorem and the correlation type reciprocity theorem. We mainly use the convolution type reciprocity theorem.

[0090] G(x|A,ω) represents the Green's function of the seismic wave response excited by point A inside the medium volume V and received by point x inside the volume V. According to formula (1), the Helmholtz equations of points A and B inside the medium can be defined, namely:

[0091]

[0092]

[0093] in, represents the gradient operator, k represents the wave number, ω represents the angular frequency, δ represents the step function, G and G0 represent the Green's function, and x, A and B represent the coordinate positions. Wave number k = ω / v(g), k0 = ω / v0(g), x∈V, since the medium is isotropic and has uniform velocity, if A and B are inside the medium, then v(x) = v0(x), and if they are outside, then v(x) ≠ v0(x). Multiply both sides of formula (2) by G0(x|B,ω), and multiply both sides of formula (3) by G(x|A,ω) and then subtract them, and we get:

[0094]

[0095] in, represents the gradient operator, ω represents the angular frequency, δ represents the step function, G and G0 represent the Green's function, and x, A and B represent the coordinate positions. Since the point x is inside the volume V, the method of integration by parts can be obtained:

[0096]

[0097]

[0098] in, represents the gradient operator, ω represents the angular frequency, G and G0 represent the Green's function, and x, A and B represent the coordinate positions.

[0099] Substituting formula (5) and formula (6) into formula (4), we can obtain:

[0100]

[0101] in, represents the gradient operator, ω represents the angular frequency, G and G0 represent the Green's function, and x, A and B represent the coordinate positions.

[0102] Using Gauss's theorem to integrate the closed surface S of volume V, we can get:

[0103]

[0104] in, represents the gradient operator, ω represents the angular frequency, G and G0 represent the Green's function, x, A and B represent the coordinate position, Represents the unit external normal vector of the boundary of surface S.

[0105] In order to simplify the convolution reciprocity theorem and the cross-correlation reciprocity theorem, so that they only contain the integral kernel composed of two monopole sources, and it is easier to understand the physical meaning of the two, the far-field approximation theory is introduced. The incident wave in the homogeneous medium will cause the scatterer to scatter, which can be expressed in the far field as:

[0106] G(x|B,ω)~m(φ,θ)e ikr / r (9)

[0107] Where ω represents the angular frequency, G represents the Green's function, x and B represent the coordinate positions, r represents the distance from B to the scattering point x in the scatterer, and m(φ,θ) is only a function of the angular coordinates and therefore has no effect on the partial derivative of r. The far field approximation of the gradient term in formula (10) is:

[0108]

[0109] in, represents the gradient operator, ω represents the angular frequency, G represents the Green's function, x and A represent the coordinate position, represents the unit external normal vector of the boundary of the surface S, r represents the distance from the scatterer A to the scattering point x, and k represents the wave number.

[0110] Because there are:

[0111]

[0112] Where ω represents the angular frequency, G represents the Green's function, x and A represent the coordinate positions, r represents the distance from A inside the scatterer to the scattering point x far away, and k represents the wave number.

[0113] So we have:

[0114]

[0115] Where ω represents the angular frequency, G represents the Green's function, x and A represent the coordinate positions, r represents the distance from A inside the scatterer to the scattering point x far away, and k represents the wave number. represents the unit external normal vector of the surface S boundary, represents the gradient operator and i represents the imaginary unit.

[0116] Similarly, there are:

[0117]

[0118] Where ω represents the angular frequency, G0 represents the Green's function, x and B represent the coordinate positions, r represents the distance from B inside the scatterer to the scattering point x far away, and k represents the wave number. represents the unit external normal vector of the surface S boundary, represents the gradient operator and i represents the imaginary unit.

[0119] Substituting formula (12) (13) into formula (8) yields:

[0120]

[0121] Among them, i represents the imaginary unit, G represents the Green's function, k represents the wave number, ω represents the angular frequency, x, A and B represent the coordinate points, S represents the integral curve, and * represents taking the conjugate.

[0122] Step 4: Use the obtained first-order multiple waves to perform conventional multiple wave imaging and multiple wave imaging:

[0123] The imaging principle of multiple wave reverse time migration is to carry out forward and reverse extension of the wave field, and then cross-correlate the two to obtain the imaging result. Only the k-1 order multiple waves and the k order multiple waves of the seismic record can be imaged by cross-correlating. The steps of multiple wave order reverse time migration are:

[0124] 1. Since the original shot record of the free surface boundary obtains a mixed wave of primary waves and multiple waves, this method uses the original shot record of the surface absorption boundary to obtain only a single wave for subsequent processing operations;

[0125] 2. Based on the primary wave record, the convolution type seismic interferometry is used to predict the first order multiple waves. On the basis of the prediction of the first order multiple waves, the second order multiple waves are predicted, and the N order multiple waves are derived from the first order;

[0126] 3. Based on the zero-order multiple waves and the first-order multiple waves, migration calculation is performed as the forward wave field and the reverse wave field to obtain the first-order multiple wave reverse time migration imaging section;

[0127] 4. Use k-1 order multiple waves and k order multiple waves as forward wave field and reverse wave field for migration calculation respectively until k = N, and the reverse time migration imaging profiles of each order of multiple waves can be obtained.

[0128] The formula for reverse time migration of multiple waves of different orders and the imaging conditions of multiple waves of different orders can be expressed as follows:

[0129]

[0130] Wherein, m represents the order of multiple waves, F represents the forward wave field, B represents the backward wave field, M represents the wave field value, x and z represent the coordinate points, and t represents the time. In this method, only the first-order multiple waves are used, so N=1.

[0131] Example 2

[0132] In the specific embodiment 2 of the present invention, a multiple wave reverse time migration technology based on the convolution interference principle of the present invention is applied to the theoretical model and achieves an ideal calculation effect.

[0133] First, test the simple model. Figure 1 A three-level horizontal model is shown. Figure 2 A three-level horizontal model single-gun record is shown, in which Figure 2 a is a primary wave; Figure 2 b is multiple waves; Figure 2 c is the full wave field. By comparison, it can be seen that this method successfully predicts multiple waves using a single wave. Figure 3 A three-layer horizontal model reverse time migration is shown, where Figure 3 a is a primary wave; Figure 3 b is the conventional multiple wave; Figure 3 b is the multiple waves of different orders. The correctness of the algorithm is verified by a simple model.

[0134] Example 3

[0135] The inclined layer model was tested. Figure 4 The oblique layer model is shown. Figure 5 The single shot record of the inclined layer model is shown, Figure 5 a is a primary wave; Figure 5 b is multiple waves; Figure 5 c is the full wave field. Figure 6 The reverse time migration of a salt dome model is shown, where Figure 6 a is a primary wave; Figure 6 b is the multiple waves of different orders. By comparison, it can be seen that the reverse time migration effect of the multiple waves obtained by this method reduces the noise in the horizontal layer.

[0136] Example 4

[0137] Then the complex salt dome model was tested. Figure 7 A salt dome model is shown. Figure 8 A single shot record of a salt dome model is shown, in which Figure 8 a is a primary wave; Figure 8 b is multiple waves; Figure 8 c is the full wave field. Fig. 9 The reverse time migration of a salt dome model is shown, where Fig. 9 a is a primary wave; Fig. 9 b is the conventional multiple wave; Fig. 9 b is the order multiple wave. By comparison, it can be seen that the reverse time migration effect of the multiple waves obtained by this method improves the imaging accuracy of the subsalt structure, as shown in the ellipse. The universality of this method is verified by a complex model.

[0138] The invention proposes a multiple wave reverse time migration technology based on the convolution interference principle, which can realize the processing of multiple wave reverse time migration in different orders, and improve the subsalt imaging accuracy of the traditional reverse time migration.

[0139] Finally, it should be noted that the above is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art can still modify the technical solutions recorded in the aforementioned embodiments or replace some of the technical features therein with equivalents. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

[0140] Except for the technical features described in the specification, all other technical features are known technologies to those skilled in the art.

Claims

1. A multiple wave reverse time migration method based on the convolution interference principle, characterized in that: The multiple wave reverse time migration method based on the convolution interference principle includes: Step 1, input velocity field and establish observation system; Step 2, using the surface absorbing boundary condition to generate a primary wave without direct waves; Step 3, performing a seismic interferometric convolution process on the primary wave to obtain a first-order multiple wave; Step 4, using the obtained first-order multiple waves to perform conventional multiple wave imaging and fractional multiple wave imaging.

2. The multiple wave reverse time migration method based on the convolution interference principle according to claim 1, characterized in that: In step 1, the source and receiver points are placed on the surface for excitation, and the shot spacing and trace spacing are set to three grid points.

3. The multiple wave reverse time migration method based on the convolution interference principle according to claim 1, characterized in that: In step 3, in 3D space, for the acoustic wave medium, the acoustic wave Green's function in the frequency domain satisfies the three-dimensional Helmholtz equation: in, represents the gradient operator, k represents the wave number, ω represents the angular frequency, δ represents the step function, s represents the source, g represents the detection point, and G represents the Green's function; Formula (1) is a second-order partial differential equation with two solutions, namely the causal Green's function G(g|s,ω) and the non-causal Green's function G(g|s,ω) * ; The positions of the excitation point and the detection point in the Green's function can be interchanged, that is, G(g|s,ω)=G(s|g,ω), which conforms to the reciprocity theorem.

4. The multiple wave reverse time migration method based on the convolution interference principle according to claim 3, characterized in that: In step 3, G(x|A,ω) represents the Green's function of the seismic wave response excited by point A inside the medium volume V and received by point x inside the volume V. According to formula (1), the Helmholtz equation of points A and B inside the medium can be defined, namely: in, represents the gradient operator, k represents the wave number, ω represents the angular frequency, δ represents the step function, G and G0 represent Green's functions, x, A and B represent the coordinate positions; wave number k = ω / v(g), k0 = ω / v0(g), x∈V, since the medium is isotropic and has uniform velocity, if A and B are inside the medium, then v(x) = v0(x), and if they are outside, then v(x)≠v0(x).

5. The multiple wave reverse time migration method based on the convolution interference principle according to claim 4, characterized in that: In step 3, multiply both sides of formula (2) by G0(x|B,ω), and multiply both sides of formula (3) by G(x|A,ω) and then subtract them to obtain: in, represents the gradient operator, ω represents the angular frequency, δ represents the step function, G and G0 represent the Green's function, x, A and B represent the coordinate position; since the point x is inside the volume V, the method of integration by parts can be obtained: in, represents the gradient operator, ω represents the angular frequency, G and G0 represent the Green's function, and x, A and B represent the coordinate positions.

6. The multiple wave reverse time migration method based on the convolution interference principle according to claim 5, characterized in that: In step 3, substitute formula (5) and formula (6) into formula (4) to obtain: in, represents the gradient operator, ω represents the angular frequency, G and G0 represent the Green's function, x, A and B represent the coordinate position; Using Gauss's theorem to integrate the closed surface S of volume V, we can get: in, represents the gradient operator, ω represents the angular frequency, G and G0 represent the Green's function, x, A and B represent the coordinate position, Represents the unit external normal vector of the boundary of surface S.

7. The multiple wave reverse time migration method based on the convolution interference principle according to claim 6, characterized in that: In step 3, the far-field approximation theory is introduced. The incident wave in the homogeneous medium will cause the scatterer to scatter, which can be expressed in the far field as: G(x|B,ω)~m(φ,θ)e ikr / r(9) Where ω represents the angular frequency, G represents the Green's function, x and B represent the coordinate positions, r represents the distance from B to the scattering point x in the scatterer, and m(φ,θ) is only a function of the angular coordinates and therefore has no effect on the partial derivative of r. The far field approximation of the gradient term in formula (10) is: in, represents the gradient operator, ω represents the angular frequency, G represents the Green's function, x and A represent the coordinate position, represents the unit external normal vector of the boundary of the surface S, r represents the distance from the scatterer A to the scattering point x, and k represents the wave number.

8. The multiple wave reverse time migration method based on the convolution interference principle according to claim 7, characterized in that: In step 3, because there are: Where ω represents the angular frequency, G represents the Green's function, x and A represent the coordinate positions, r represents the distance from A inside the scatterer to the scattering point x far away, and k represents the wave number; So we have: Where ω represents the angular frequency, G represents the Green's function, x and A represent the coordinate positions, r represents the distance from A inside the scatterer to the scattering point x far away, and k represents the wave number. represents the unit external normal vector of the surface S boundary, represents the gradient operator, i represents the imaginary unit; Similarly, there are: Where ω represents the angular frequency, G0 represents the Green's function, x and B represent the coordinate positions, r represents the distance from B inside the scatterer to the scattering point x far away, and k represents the wave number. represents the unit external normal vector of the surface S boundary, represents the gradient operator and i represents the imaginary unit.

9. The multiple wave reverse time migration method based on the convolution interference principle according to claim 8, characterized in that: In step 3, substitute formula (12) (13) into formula (8) to obtain: Among them, i represents the imaginary unit, G represents the Green's function, k represents the wave number, ω represents the angular frequency, x, A and B represent the coordinate points, S represents the integral curve, and * represents taking the conjugate.

10. The multiple wave reverse time migration method based on the convolution interference principle according to claim 1, characterized in that: In step 4, the steps of multiple-order RTM are: Step 41. Since the original shot record of the free surface boundary obtains a mixed wave of primary waves and multiple waves, the original shot record of the surface absorption boundary is used to obtain only a primary wave for subsequent processing operations; Step 42. Based on the primary wave record, apply convolution type seismic interferometry to predict the first order multiple waves, predict the second order multiple waves based on the predicted first order multiple waves, and deduce the N order multiple waves from the first order; Step 43. Based on the zero-order multiple waves and the first-order multiple waves, migration calculation is performed as the forward wave field and the reverse wave field to obtain the first-order multiple wave reverse time migration imaging section; Step 44. Use k-1 order multiple waves and k order multiple waves as forward wave field and reverse wave field to perform migration calculation respectively until k=N, and the reverse time migration imaging sections of each order of multiple waves can be obtained.

11. The multiple wave reverse time migration method based on the convolution interference principle according to claim 10, characterized in that: In step 4, the formula for reverse time migration of multiple waves of different orders and the imaging conditions of multiple waves of different orders can be expressed as: Wherein, m represents the order of multiple waves, F represents the forward wave field, B represents the reverse wave field, M represents the wave field value, x and z represent the coordinate points, and t represents the time. Here, only the first-order multiple waves are used, so N=1.

Images