Seismic reverse time migration method and system for multi-axis wave field separation based on analytic wave field

Abstract

The invention provides a seismic reverse time migration method and system based on multi-axis wave field separation of an analytic wave field, and belongs to the technical field of seismic imaging in oil-gas exploration and development. According to the method, a constructed analysis wave field is utilized to realize decomposition of the wave field, then an imaging condition is divided into a plurality of subitems, and then a final image is generated by utilizing the imaging condition. According to the method, the multi-axis decomposition of the wave field is efficiently and accurately realized based on the construction of the analytic wave field, the traditional cross-correlation imaging condition is divided into a plurality of sub-images, the offset illusion can be effectively identified and removed, the horizontal layer, the large-dip-angle layer and the fault are respectively imaged, the sub-images can be respectively processed and explained, and the imaging accuracy is improved. And the accuracy of the image and the flexibility of imaging processing are improved. In addition, a weight factor is introduced into an imaging condition, so that the resolution of a fine structure is improved, and the quality of an overall offset image is improved.

Description

Technical Field

[0001] The present invention belongs to the technical field of seismic imaging in oil and gas exploration and development, and in particular relates to a seismic reverse time migration method and system based on multi-axis wave field separation of analytical wave fields. Background Art

[0002] After investigating a large number of detectable documents and patents at home and abroad, and fully understanding and studying the software products of research institutions and geophysical companies, the existing technologies related to the present invention are summarized as follows:

[0003] Traditional reverse time migration uses a two-way wave equation to cross-correlate the source wave field obtained by the forward seismic wavelet and the receiving wave field obtained by the reverse seismic data to obtain a migration image. It can fully consider the lateral change of velocity and avoid the inclination angle limitation of the traditional migration method, and has been widely used in the industry. However, in the cross-correlation process, due to the bidirectional propagation of the wave field, the unrelated wave field components will obtain the migration artifact, which seriously affects the resolution and reliability of the imaging. Using Laplace filtering is an effective method, but it is a non-physical means that will change the phase amplitude characteristics of the image, and the effect of removing artifacts is also limited. The commonly used method at present is to separate the uplink and downlink wave fields, and eliminate the uplink source wave field and the downlink receiving wave field.

[0004] The commonly used separation method is to use the wave field vector (Poynting vector) for wave field separation, but it is based on the product of time and space derivatives. The separated wave field energy is discontinuous and the decomposition is inaccurate, which will introduce new offset errors. Decomposition in the frequency-wavenumber domain can obtain more accurate and continuous separation results, which is based on the sign of the apparent propagation velocity. However, it is very difficult to calculate the frequency of the wave field in the time domain, which requires saving the wave field values ​​of all time steps. Constructing an analytical wave field can effectively solve this problem. It constructs a complex variable analytical wave field with non-negative frequency, and wave field separation can be achieved only in the wavenumber domain.

[0005] The current wavefield separation is limited to the depth direction, which can remove some migration artifacts. However, there will still be residual artifacts in the image, which come from the non-causal horizontal traveling wave field cross-correlation, and such artifacts cannot be effectively removed by uniaxial (depth) decomposition alone. In addition, in conventional migration, structures such as horizontal planes, high-angle planes, and faults are uniformly imaged in the same result, which also increases the difficulty of further image processing and reduces the reliability of the processing and interpretation results of certain structures. Summary of the invention

[0006] The purpose of the present invention is to solve the above-mentioned problems existing in the prior art, and to provide a seismic reverse time migration method and system based on multi-axis wave field separation of analytical wave fields, so as to solve the problems of many false images, low image signal-to-noise ratio and resolution in traditional seismic reverse time migration. The method significantly improves the flexibility of the reverse time migration image processing process, the resolution and reliability of the migration image, and improves the imaging quality. The decomposed image obtained can also provide a reliable reference for subsequent migration velocity analysis and reservoir description work.

[0007] The present invention is achieved through the following technical solutions:

[0008] The first aspect of the present invention provides a seismic reverse time migration method based on multi-axis wavefield separation of analytical wavefield, wherein the method utilizes the constructed analytical wavefield to realize wavefield decomposition, and then divides the imaging conditions into multiple sub-items, and then generates the final image using the imaging conditions.

[0009] Preferably, the method comprises:

[0010] Step 1: Input seismic data and seismic wavelets;

[0011] Step 2: Construct analytical source wavelet and analytical source data;

[0012] Step 3: Obtain the source wavefield and receiving wavefield using analytical source wavelet and analytical source data;

[0013] Step 4: Obtain traveling wave fields in different directions;

[0014] Step 5: Obtaining the traveling wave field for imaging;

[0015] Step 6: Introduce weight factors into the imaging conditions to generate the final image.

[0016] Preferably, the operation of step 2 includes:

[0017] The following formula is used to perform Hilbert transform on seismic data and seismic wavelets to obtain analytical source wavelets and analytical source data:

[0018]

[0019]

[0020] in, To analyze the source wavelet, To analyze the earthquake source data, f(t) is the earthquake wavelet, d(t, x) is the earthquake data, is the Hilbert transform of the time variable.

[0021] Preferably, the operation of step three includes:

[0022] The source wave field and receiving wave field are obtained using the following formula:

[0023]

[0024]

[0025] in, represent the source wave field and the receiving wave field respectively, c is the acoustic wave velocity, δ(x-x0) represents the Dirac-Delta function, x represents the spatial position, x0 is the source position, t represents the time, ▽ 2 represents the Laplace operator.

[0026] Preferably, the operation of step 4 includes:

[0027] The source wave field and the receiving wave field are decomposed using the following formulas to obtain 8 traveling wave fields:

[0028]

[0029]

[0030]

[0031]

[0032] Where P is the acoustic wave field in the wave number domain, is the acoustic wave analytical wave field in the wave number domain, (k x ,k z ) are the components of the two-dimensional wave number, and u, d, r, and l represent up, down, left, and right, respectively.

[0033] Preferably, the operation of step five includes:

[0034] The lower left-traveling wave field of the source wave field in the wave number domain, the lower right-traveling wave field of the source wave field in the wave number domain, the upper left-traveling wave field of the receiving wave field in the wave number domain, and the upper right-traveling wave field of the receiving wave field in the wave number domain are selected from the 8 traveling waves, and then the lower left-traveling wave field of the source wave field is obtained by inverse Fourier transform. Downward right wave field of earthquake source wave field Receive the upper left wave field of the wave field Receive the upper right wave field of the wave field

[0035] Preferably, the operation of step six includes:

[0036] The final image I(x) is obtained using the following formula:

[0037]

[0038] Among them, W i represents the weight factor of a single item,

[0039] A second aspect of the present invention provides a seismic reverse time migration system based on multi-axis wavefield separation of analytical wavefields, the system comprising:

[0040] Input unit, used to input seismic data and seismic wavelets:

[0041] A construction unit, connected to the input unit, for constructing analytical source wavelets and analytical source data;

[0042] A wave field acquisition unit, connected to the construction unit, is used to obtain a source wave field and a receiving wave field by using the analytical source wavelet and the analytical source data;

[0043] A traveling wave acquisition unit, connected to the wave field acquisition unit, is used to obtain traveling wave fields in different directions;

[0044] An imaging wave field acquisition unit, connected to the traveling wave acquisition unit, for obtaining a traveling wave field for imaging;

[0045] The imaging unit is connected to the imaging wave field acquisition unit and is used to introduce weight factors into the imaging conditions to generate a final image.

[0046] Preferably, the imaging unit performs the following operations:

[0047] The final image I(x) is obtained using the following formula:

[0048]

[0049] Among them, W i represents the weight factor of a single item, represents the downward left-traveling wave field of the source wave field, represents the right-going wave field of the source wave field, Indicates the upper left wave field of the received wave field, Indicates the upper right-going wave field of the received wave field.

[0050] According to a third aspect of the present invention, a computer-readable storage medium is provided, wherein the computer-readable storage medium stores at least one computer-executable program, and when the at least one program is executed by the computer, the computer executes the steps of the seismic reverse time migration method based on multi-axis wave field separation of analytical wave fields according to the present invention.

[0051] Compared with the prior art, the beneficial effects of the present invention are as follows: based on the construction of analytical wave fields, the present invention efficiently and accurately realizes the multi-axis decomposition of wave fields, and by dividing the traditional cross-correlation imaging conditions into multiple sub-images, it can effectively identify and remove the offset artifacts, and for the horizontal planes and high-angle planes and faults, the sub-images can be processed and interpreted separately, which is conducive to improving the accuracy of the image and the flexibility of imaging processing. In addition, the present invention introduces weight factors in the imaging conditions, which is conducive to improving the resolution of fine structures, thereby improving the quality of the overall offset image. The present invention can perform reverse time offset imaging of complex underground structures, and is suitable for high-precision imaging of complex geological structures such as faults, high-angle planes, and abnormal bodies, providing strong support for subsequent offset velocity analysis and reservoir description work. BRIEF DESCRIPTION OF THE DRAWINGS

[0052] Figure 1a Lake sub-wave field;

[0053] Figure 1b After Hilbert transform;

[0054] Figure 1c Analyze wave fields;

[0055] Figure 2a Sound wave velocity model;

[0056] Figure 2b Traditional RTM results;

[0057] Figure 2c Traditional RTM plus Laplace filtering results;

[0058] Figure 3 Multi-axis traveling wave decomposition sub-images;

[0059] Figure 4a Traveling wave decomposition RTM direct stacking profile;

[0060] Figure 4b Traveling wave decomposition RTM weighted stacking profile;

[0061] Figure 5 A flowchart of the steps of the method of the present invention. DETAILED DESCRIPTION

[0062] The present invention is further described in detail below in conjunction with the accompanying drawings:

[0063] The present invention provides a seismic reverse time migration method based on multi-axis wave field separation of analytical wave fields, which is mainly used to improve the resolution of underground fine structures and remove migration artifacts to improve the flexibility of reverse time migration imaging processing and the resolution and reliability of the results.

[0064] The derivation process of the equation used in the present invention is as follows:

[0065] Step 1: Construct the wave equation for analyzing the wave field:

[0066] The commonly used acoustic wave equation can be expressed as:

[0067]

[0068] Where p is the acoustic wave field, c is the acoustic wave velocity, f(t) is the seismic wavelet, x0 is the source location, t represents time, ▽ 2 represents the Laplace operator, and δ(x-x0) represents the Dirac-Delta function.

[0069] Perform Hilbert transform on both ends of formula (1) to obtain

[0070]

[0071] in, is the Hilbert transform of the time variable. Using the linear property of the Hilbert transform and considering the independence of time and space variables, formula (2) can be expressed as:

[0072]

[0073] Combining formula (1) and formula (3), we can get the following wave equation for analytical acoustic wave field:

[0074]

[0075] in To analyze the wave field, To resolve the source wavelet. One of the advantages of the analytical wavefield is that the frequency is non-negative. Figure 1a is a 20Hz Ricker wavelet pulse wave field, whose spectrum is symmetrical about the origin. The Hilbert transform only changes the phase characteristics but not the frequency characteristics (such as Figure 1b As shown), after constructing the analytical wave field (as shown Figure 1c As shown in the figure, the amplitude of the negative part of the spectrum is 0, which greatly facilitates the subsequent wave field separation.

[0076] It should be noted that in reverse time migration, if the source wave field adopts an analytical wave field, the return data also needs to be Hilbert transformed accordingly to meet the corresponding analytical conditions. That is, when generating the receiving wave field, formula (4) is also used, except that the right side Represents the parsed return data. The method of generating the parsed return data is to use the original data as the real part and its Hilbert transform as the imaginary part.

[0077] Step 2: Construct the traveling wave decomposition wave field. Traveling wave decomposition is usually based on the sign of the apparent propagation velocity, which is defined as v j =ω / kj (where v j is the apparent propagation velocity, ω is the angular frequency, k j is the wave number component), where j = x, y, z (representing three coordinate axes respectively). Since the calculation frequency needs to store the wave field of each time step during the wave field propagation in the time-space domain, the calculation efficiency is extremely low. The above analytical wave field can avoid this problem. Using the complex variable analytical wave field, the two-dimensional traveling wave field can be divided into four parts along the x and z axes:

[0078]

[0079]

[0080]

[0081]

[0082] Where P is the acoustic wave field in the wave number domain, is the acoustic wave analytical wave field in the wave number domain, (k x ,k z ) are the components of the two-dimensional wave number, and u, d, r, and l represent the up, down, left, and right traveling waves, for example In other words, the three-dimensional wave field can be divided into eight parts according to the direction of the component. For the convenience of description and calculation, the following two-dimensional wave field is used as an example. is the wave field in the time-wavenumber domain, It is the wave field in the time-space domain.

[0083] Step 3: Introduce weight factors into the imaging conditions. The traditional cross-correlation imaging conditions can be written as follows:

[0084]

[0085] Where S = (S d ,S u ) T and R = (R d ,R u ) T are the source wave field and the receiving wave field. For the multi-axis decomposed wave field, formula (6) becomes the sum of 16 terms, and a weight matrix is ​​introduced, so the imaging condition can be written as:

[0086]

[0087] Among them, W is the weight matrix.

[0088] In formula (6), S d R d , Su R u and S u R d Usually an offset artifact, S d R u This is usually the only image retained, so the weight matrix is ​​modified to

[0089]

[0090] Where W i Represents the weight of a single item (1, 2, 3, 4 represent four coefficients, which have no specific physical meaning), and their sum is 1. Therefore, formula (7) can be simplified as:

[0091]

[0092] in and are the same direction related items, and is the reverse correlation term. In reverse time migration, the same-direction correlation term can often characterize flat layers, while the reverse correlation term can often characterize inclined layers and complex structures, and some illusions will remain. Formula (9) only selects the downward source wavefield and the upward receiving wavefield. Therefore, the present invention removes the non-causal wavefield cross-correlation of the upward source wavefield and the downward receiving wavefield through traveling wave decomposition, removes the migration illusion from the perspective of wavefield propagation, and can further remove the migration illusion for the sub-image through directional filtering.

[0093] The present invention utilizes the constructed analytical wave field to achieve accurate decomposition of the wave field. Then, by dividing the traditional imaging conditions into multiple sub-items, the horizontal plane and the high-angle plane and the fault are imaged separately, and the obtained sub-items can be imaged separately, thereby improving the flexibility of imaging.

[0094] The present invention introduces weight factors into traditional imaging conditions, effectively removes offset artifacts and improves the resolution of fine structures, thereby improving the overall imaging quality and reliability, and providing a reference for subsequent processing and analysis.

[0095] Specifically, Figure 5 As shown, the method of the present invention comprises:

[0096] Step 1: Input earthquake data and earthquake wavelets:

[0097] Step 2: Construct analytical source wavelet and analytical source data;

[0098] Step 3: Obtain the source wavefield and receiving wavefield using analytical source wavelet and analytical source data;

[0099] Step 4: Obtain traveling wave fields in different directions;

[0100] Step 5: Obtaining the traveling wave field for imaging;

[0101] Step 6: Introduce weight factors into the imaging conditions to generate the final image.

[0102] The embodiments of the method of the present invention are as follows:

[0103] Embodiment 1:

[0104] The operation of step one includes: obtaining seismic data and seismic wavelets from conventional seismic data.

[0105] Embodiment 2:

[0106] The operations of step 2 include:

[0107] The following formula is used to perform Hilbert transform on seismic data and seismic wavelets to obtain analytical source wavelets and analytical source data:

[0108]

[0109]

[0110] in To analyze the source wavelet, To analyze the earthquake source data, f(t) is the earthquake wavelet, d(t, x) is the earthquake data, is the Hilbert transform of the time variable.

[0111] Embodiment 3:

[0112] The operations in step three include:

[0113] The source wave field and receiving wave field are obtained using the following formula:

[0114]

[0115]

[0116] in, represent the source wave field and the receiving wave field respectively, c is the acoustic wave velocity, δ(x-x0) represents the Dirac-Delta function, x represents the spatial position, x0 is the source position, t represents the time, ▽ 2 represents the Laplace operator.

[0117] Embodiment 4:

[0118] The operations in step 4 include:

[0119] The source wave field and the receiving wave field are decomposed using the following formulas to obtain 8 traveling wave fields:

[0120]

[0121]

[0122]

[0123]

[0124] Where P is the acoustic wave field in the wave number domain, is the acoustic wave analytical wave field in the wave number domain, (k x ,k z ) are the components of the two-dimensional wave number, u, d, r, and l represent the upward, downward, left, and right traveling waves, respectively. Can be written as They are spatially Fourier transformed to obtain the formula (5): The traveling wave decomposition procedure is consistent for both the source wavefield and the received wavefield.

[0125] Embodiment 5:

[0126] The operations of step five include:

[0127] The lower left-traveling wave field of the source wave field in the wave number domain, the lower right-traveling wave field of the source wave field in the wave number domain, the upper left-traveling wave field of the receiving wave field in the wave number domain, and the upper right-traveling wave field of the receiving wave field in the wave number domain are selected from the 8 traveling waves, and then the lower left-traveling wave field of the source wave field is obtained by inverse Fourier transform. Downward right wave field of earthquake source wave field Receive the upper left wave field of the wave field Receive the upper right wave field of the wave field

[0128] and These four parameters correspond to the inverse Fourier transform of the wave field shown in formula (5), for example, The left-going wave field representing the source wave field is obtained by decomposing the source wave field into After inverse Fourier transform, we get The upper right wave field representing the received wave field is obtained by decomposing the received wave field into Obtained by inverse Fourier transform.

[0129] Embodiment 6:

[0130] The operations of step six include:

[0131] The final image I(x) is obtained using the following formula:

[0132]

[0133] From formula (9), we can see that the present invention only uses and Perform cross-correlation.

[0134] In formula (9), W is a weight factor, and its value can be adjusted according to the actual effect. The present invention adopts the structure similarity index of the fitting formula (9) and the traditional cross-correlation imaging condition, and uses the gradient method to invert and obtain the value of the weight factor W.

[0135] in, and are the same direction related items, and is the reverse correlation term. i Represents the weight factor of a single item, and the sum of all weight factors is 1, that is,

[0136] The embodiments of the present invention are as follows:

[0137] Embodiment 7:

[0138] A typical Marmousi model was selected for migration imaging. The model area is 7.5 km × 3 km, the spatial spacing is 10 m, and the seismic wavelet is a 20 Hz Ricker wavelet. Figure 2a to Figure 2b Results from conventional reverse time migration imaging are shown. Figure 2a The speed model. Figure 2b For traditional RTM results, the migration artifacts are very serious and the imaging quality of underground structures is very poor. After using Laplace filtering, the low-frequency migration artifact problem has been alleviated to a certain extent, but some high-frequency noise still exists.

[0139] Figure 3 In the sub-images of the traveling wave decomposition, the same-direction correlation term can characterize the flat layer, while the reverse correlation term can characterize the large dip angle and fault. At the same time, some false images will be retained, but these false images are in the flat layer area and can be well identified and removed. Here, a directional filter can be used to remove them.

[0140] Figure 4a and Figure 4b is the reconstructed image. Figure 4a It is a simple superposition of four sub-images. Although the offset artifacts are reduced, the amplitudes of the large-angle planes and faults are small. Figure 4b It is the weighted superposition of the present invention, and the weights are obtained by inversion, (W1, W2, W3, W4) = (0.160, 0.293, 0.387, 0.160). It can be seen Figure 4b The faults and structures are clearer and the resolution is higher.

[0141] The present invention also provides a seismic reverse time migration system based on multi-axis wave field separation of analytical wave fields, and the embodiments of the system are as follows:

[0142] Embodiment 8:

[0143] The system comprises:

[0144] Input unit, used to input seismic data and seismic wavelets:

[0145] A construction unit, connected to the input unit, for constructing analytical source wavelets and analytical source data;

[0146] A wave field acquisition unit, connected to the construction unit, is used to obtain a source wave field and a receiving wave field by using the analytical source wavelet and the analytical source data;

[0147] A traveling wave acquisition unit, connected to the wave field acquisition unit, is used to obtain traveling wave fields in different directions;

[0148] An imaging wave field acquisition unit, connected to the traveling wave acquisition unit, for obtaining a traveling wave field for imaging;

[0149] The imaging unit is connected to the imaging wave field acquisition unit and is used to introduce weight factors into the imaging conditions to generate a final image.

[0150] Specifically, the imaging unit performs the following operations:

[0151] The final image I(x) is obtained using the following formula:

[0152]

[0153] Among them, W i represents the weight factor of a single item, represents the downward left-traveling wave field of the source wave field, represents the right-going wave field of the source wave field, Indicates the upper left wave field of the received wave field, Indicates the upper right-going wave field of the received wave field.

[0154] The present invention also provides a computer-readable storage medium, which stores at least one computer-executable program. When the at least one program is executed by the computer, the computer executes the steps of the seismic reverse time migration method based on multi-axis wave field separation of analytical wave fields according to the present invention.

[0155] The present invention is mainly used to improve the reliability and resolution of seismic reverse time migration imaging, and can provide strong support for subsequent migration velocity analysis and reservoir description work.

[0156] The present invention has a broad application prospect in the exploration and development of oil and gas fields. First, in terms of technology, the high-precision imaging of reverse time migration based on analytical wave fields is still being improved and developed, and the present invention can provide theoretical support for the research of related algorithms; secondly, in terms of application, the present invention can be applied in actual production, and has a relatively good application prospect in complex exploration areas such as Northwest and North China of Sinopec, and can provide effective technical support for the migration velocity analysis and reservoir description work in related exploration areas.

[0157] In the description of the present invention, it should be noted that, unless otherwise clearly specified and limited, the terms "connected" and "connection" should be understood in a broad sense, for example, it can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection or an electrical connection; it can be a direct connection or an indirect connection through an intermediate medium. For ordinary technicians in this field, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.

[0158] In the description of the present invention, unless otherwise specified, the terms "upper", "lower", "left", "right", "inside", "outside", etc. indicate directions or positional relationships based on the directions or positional relationships shown in the accompanying drawings. They are only for the convenience of describing the present invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific direction, be constructed and operated in a specific direction. Therefore, they cannot be understood as limitations on the present invention.

[0159] The above technical solution is only one implementation mode of the present invention. For those skilled in the art, it is easy to make various types of improvements or modifications based on the principles disclosed in the present invention, and it is not limited to the technical solution described in the above specific embodiments of the present invention. Therefore, the above description is only preferred and does not have a restrictive meaning.

Claims

1. A seismic reverse time migration method based on multi-axis wave field separation of analytical wave fields, characterized by: The method utilizes the constructed analytical wave field to realize the decomposition of the wave field, then divides the imaging condition into a plurality of sub-items, and then utilizes the imaging condition to generate the final image.

2. The seismic reverse time migration method based on multi-axis wavefield separation of analytical wavefield according to claim 1, characterized in that: The method comprises: Step 1: Input seismic data and seismic wavelets; Step 2: Construct analytical source wavelet and analytical source data; Step 3: Obtain the source wavefield and receiving wavefield using analytical source wavelet and analytical source data; Step 4: Obtain traveling wave fields in different directions; Step 5: Obtaining the traveling wave field for imaging; Step 6: Introduce weight factors into the imaging conditions to generate the final image.

3. The seismic reverse time migration method based on multi-axis wavefield separation of analytical wavefield according to claim 2, characterized in that: The operations of step 2 include: The following formula is used to perform Hilbert transform on seismic data and seismic wavelets to obtain analytical source wavelets and analytical source data: in, To analyze the source wavelet, To analyze the earthquake source data, f(t) is the earthquake wavelet, d(t, x) is the earthquake data, is the Hilbert transform of the time variable.

4. The seismic reverse time migration method based on multi-axis wavefield separation of analytical wavefield according to claim 3, characterized in that: The operations in step three include: The source wave field and receiving wave field are obtained using the following formula: in, represent the source wave field and the receiving wave field respectively, c is the acoustic wave velocity, δ(x-x0) represents the Dirac-Delta function, x represents the spatial position, x0 is the source position, t represents the time, ▽ 2 represents the Laplace operator.

5. The seismic reverse time migration method based on multi-axis wave field separation of analytical wave field according to claim 4, characterized in that: The operations in step 4 include: The source wave field and the receiving wave field are decomposed using the following formulas to obtain 8 traveling wave fields: Where P is the acoustic wave field in the wave number domain, is the acoustic wave analytical wave field in the wave number domain, (k x ,k z ) are the components of the two-dimensional wave number, and u, d, r, and l represent up, down, left, and right, respectively.

6. The seismic reverse time migration method based on multi-axis wave field separation of analytical wave fields according to claim 5, characterized in that: The operations of step five include: The lower left-traveling wave field of the source wave field in the wave number domain, the lower right-traveling wave field of the source wave field in the wave number domain, the upper left-traveling wave field of the receiving wave field in the wave number domain, and the upper right-traveling wave field of the receiving wave field in the wave number domain are selected from the 8 traveling waves, and then the lower left-traveling wave field of the source wave field is obtained by inverse Fourier transform. Downward right wave field of earthquake source wave field Receive the upper left wave field of the wave field Receive the upper right wave field of the wave field 7. The seismic reverse time migration method based on multi-axis wavefield separation of analytical wavefield according to claim 6, characterized in that: The operations of step six include: The final image I(x) is obtained using the following formula: Among them, W i represents the weight factor of a single item, 8. A seismic reverse time migration system based on multi-axis wave field separation of analytical wave fields, characterized by: The system comprises: Input unit, used to input seismic data and seismic wavelets: A construction unit, connected to the input unit, for constructing analytical source wavelets and analytical source data; A wave field acquisition unit, connected to the construction unit, is used to obtain a source wave field and a receiving wave field by using an analytical source wavelet and analytical source data; A traveling wave acquisition unit, connected to the wave field acquisition unit, is used to obtain traveling wave fields in different directions; An imaging wave field acquisition unit, connected to the traveling wave acquisition unit, for obtaining a traveling wave field for imaging; The imaging unit is connected to the imaging wave field acquisition unit and is used to introduce weight factors into the imaging conditions to generate a final image.

9. The seismic reverse time migration system based on multi-axis wavefield separation of analytical wavefield according to claim 8, characterized in that: The imaging unit performs the following operations: The final image I(x) is obtained using the following formula: Among them, W i represents the weight factor of a single item, represents the left-traveling wave field of the source wave field, represents the right-going wave field of the source wave field, Indicates the upper left wave field of the received wave field, Indicates the upper right-going wave field of the received wave field.

10. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores at least one computer-executable program, and when the at least one program is executed by the computer, the computer executes the steps of the seismic reverse time migration method based on multi-axis wavefield separation of analytical wavefields according to any one of claims 1 to 7.

Images