A phase-controllable and efficient reverse time migration imaging method and device

Abstract

The application provides a phase-controllable high-efficiency reverse-time migration imaging method, and mainly solves the problem that the existing reverse-time imaging method cannot solve the wavelet space variation and complex wave field in post-stack phase shaping processing. The application comprises the following steps: S1. obtaining shot gather data after fidelity processing and a high-precision depth domain velocity model; S2. adjusting the phase of a source wave field input seismic wavelet; S3. for any shot seismic data, completing optimized time source wave field forward continuation and total length time geophone point wave field reverse-time continuation and applying an imaging condition, completing reverse-time migration processing of any shot data, and obtaining reverse-time imaging data volume and amplitude energy data volume of the shot; and S4. obtaining reverse-time migration data volume and amplitude energy data volume of the whole work area. According to the type of the seismic wave equation adopted, the corresponding source wavelet phase is changed, the phase of the imaging processing result is adaptively changed in the reverse-time migration process, and the accurate and fine description of a complex reservoir or a complex structure is facilitated.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic exploration data processing, and in particular to a phase-controllable and efficient reverse time migration imaging method and device. Background Art

[0002] Seismic wave time-reversal imaging (RTI) uses fewer approximations to the seismic wave equation, can resolve multivalued traveltime problems, and is suitable for imaging at arbitrarily steep inclination angles and in situations with dramatic lateral and vertical velocity variations. It can also accurately image wave fields, such as multiples and gyratory reflections, which are often considered interference waves. Compared to Kirchhoff prestack depth imaging and single-pass seismic depth imaging, RTI offers significant advantages, making it one of the most theoretically mature and accurate high-precision seismic imaging technologies currently available.

[0003] However, a problem with reverse time migration imaging is that the current reverse time migration imaging results have significant differences in wave group phase characteristics compared to conventional Kirchhoff migration results. Some wave group characteristics are similar, while others even have phase differences of 90 degrees. If the data after migration stacking is subjected to post-stack phase shaping, the non-stretching distortion caused by the migration and the mutual interference of complex wave fields in complex structures are already superimposed on the post-stack imaging data. Therefore, it is no longer possible to address the differences in wave group characteristics introduced by the sub-wave space variation and observation orientation during the migration process on the post-stack data, which is not conducive to the accurate and detailed characterization of complex reservoirs or complex structures. At the same time, the existing reverse time migration process uses the entire time length of the seismic record for reverse time imaging, in which some wave field imaging has little effect on the accuracy of the final imaging result, so there is still much room for improvement in imaging efficiency. Summary of the Invention

[0004] The technical problem to be solved by the present invention is to overcome the problem of prior art that imaging stacking results cannot correct the phase of deep and space-time variations of deep-domain seismic wavelets. Instead, a phase-controllable and efficient reverse time migration imaging method is provided. This phase-controllable and efficient reverse time migration imaging method can adapt the phase of the corresponding source wavelet according to the type of seismic wave equation used, effectively achieving adaptive phase changes in the imaging processing results during the reverse time migration process. This effectively solves the problems of wavelet space-time variation and wavefield complexity that exist in post-stack phase shaping processing, thereby facilitating the accurate and detailed characterization of complex reservoirs or complex structures. The present invention also provides a phase-controllable reverse time migration imaging device.

[0005] The present invention solves the problem by the following technical solution: The phase-controllable and efficient reverse time migration imaging method comprises the following steps:

[0006] S1. Obtain shot gather data and high-precision depth-domain velocity models after fidelity processing in the seismic area;

[0007] S2. adjusting the phase of the input seismic wavelet that generates the source wavefield according to the type of seismic wave equation used;

[0008] S3. For any shot seismic data in the shot data in step S1, the forward extension of the source wave field of the optimized time and the reverse extension of the total track long-time detection point wave field are completed;

[0009] S4. At any time step optimized in step S3, applying the imaging condition to any imaging point in the imaging space of the shot data, and accumulating the seismic wavefield after the imaging condition is applied within the 0-T time, thereby completing the reverse time migration imaging processing of the shot data in the shot gather data;

[0010] S5. Complete the reverse time migration imaging processing of any single shot data in the shot gather data in step S1 according to step S3 and step S4 to obtain the single shot reverse time imaging data volume and amplitude energy data volume;

[0011] S6. Complete the reverse time migration processing of all shot data in the work area according to step S5, and obtain the reverse time migration data volume and amplitude energy data volume of the entire work area by cutting and superimposing all single-shot reverse time imaging data volumes and amplitude energy data volumes, and complete the reverse time migration optimization processing.

[0012] Preferably, the shot gather data in step S1 needs to be subjected to surface consistent deconvolution processing and surface consistent amplitude energy compensation processing to ensure that the seismic wavelet amplitude energy, phase, and frequency of all shot points and receiver points are spatially consistent;

[0013] The surface consistent deconvolution processing and surface consistent energy compensation processing both use a relatively long time window, at least including the target layer segment for more than 1 second, to ensure that the processed seismic data can retain the complete wave field propagation effects and energy relationships such as seismic reflection and transmission.

[0014] Preferably, the method for obtaining the high-precision depth-domain velocity model is obtained by adopting a sophisticated three-dimensional velocity inversion method, specifically including a full waveform inversion method and a multi-information constrained grid tomography velocity inversion method;

[0015] The multi-information constrained grid tomographic velocity inversion method obtains a velocity model representing fine structures by iterating the depth domain velocity model, and specifically comprises the following steps:

[0016] 11) Use prestack time migration velocity to perform time-to-depth conversion, and use VSP velocity or logging velocity, interpreted horizon information, and other information to constrain the construction of an initial depth domain velocity model. Complete the first round of prestack depth migration, output the first round of depth domain imaging volume and depth domain imaging gathers, and use a migration grid of 10 bins or less for the main and tie lines.

[0017] 12) Completing the residual delay spectrum calculation and seed point picking for the first round of depth domain imaging gathers, with the depth picking density being less than or equal to 40 times the depth step; and completing the three-dimensional grid tomographic velocity inversion, with the spatial inversion grid being less than or equal to 20 bins and the depth inversion grid being less than or equal to 40 times the depth step. The three-dimensional grid tomographic inversion velocity model in this step is added to the initial depth domain model in step 11) to obtain a velocity model representing large-scale structural features;

[0018] 13) Complete the second round of pre-stack depth migration and output the second round of depth domain imaging volume and depth domain imaging gathers. The migration grid of the main survey line and the tie survey line should be less than or equal to 5 bins.

[0019] 14) Completing the residual delay spectrum calculation and seed point picking for the second round of depth domain imaging gathers, with the picking density in the depth direction being less than or equal to 20 times the depth step; and completing the three-dimensional grid tomographic velocity inversion, with the spatial inversion grid being less than or equal to 10 bins and the depth inversion grid being less than or equal to 20 times the depth step. The three-dimensional grid tomographic inversion velocity model in this step is added to the velocity model representing large-scale structural features in step 12), thereby obtaining a velocity model representing medium- and large-scale complex structures;

[0020] 15) Complete the third round of pre-stack depth migration and output the third round of depth domain imaging volume and depth domain imaging gathers. The migration grid of the main survey line and the tie survey line is less than or equal to 2 bins.

[0021] 16) Completing the residual delay spectrum calculation and seed point picking for the third round of depth domain imaging gathers, with the picking density in the depth direction being less than or equal to 20 times the depth step; and completing the three-dimensional grid tomographic velocity inversion, with the spatial inversion grid being less than or equal to 10 bins and the depth inversion grid being less than or equal to 20 times the depth step. The three-dimensional grid tomographic inversion velocity model in this step is added to the velocity model representing the large-scale complex structures in step 14) to obtain a velocity model representing the fine structure;

[0022] 17) When the depth-domain imaging data volume corresponding to the depth-domain velocity model meets the geological requirements, the velocity model matches the imaging result structure, and the imaging gather flattening standard, a velocity model that characterizes the fine structure that meets the standard is obtained; and the iteration of the depth-domain velocity model ends;

[0023] Preferably, the factors that satisfy the geological requirements in step 17) are accurate structural imaging position, clear geological body imaging, and clear fault and fracture depiction;

[0024] The flattening standard requirement of the imaging gather is that the energy cluster on the residual delay spectrum of the depth domain imaging gather is returned to zero.

[0025] The conditions for the structure matching between the velocity model and the imaging results are: the top and bottom depth positions of the velocity changes in the depth domain velocity model are consistent with the depth position and trend of the marker layer in the seismic imaging results; the local velocity changes in the depth domain velocity model are consistent with the local geological body change range in the seismic imaging results.

[0026] Preferably, the seismic wavelet for generating the source wavefield in step S2 includes a zero-phase Ricker wavelet or a zero-phase frequency-domain band-limited seismic wavelet;

[0027] The method for adjusting the phase of the input seismic wavelet for generating the source wavefield using the seismic wave equation type adopted in step S2 comprises the following steps:

[0028] 21) To make the wave group characteristics of the reverse time imaging results consistent with those of the Kirchhoff depth imaging results, it is necessary to perform a 90-degree phase rotation on the input seismic wavelet that generates the source wavefield;

[0029] 22) If the reverse time imaging result is to have a wave group characteristic close to zero phase, the input seismic wavelet generating the source wave field is not subjected to phase rotation;

[0030] 23) To make the reverse time imaging result have The wave group characteristics of the phase can be obtained by Phase rotation is achieved;

[0031] Preferably, the types of seismic wave equations used include first-order seismic wave equations and second-order seismic wave equations;

[0032] Taking the three-dimensional acoustic wave equation as an example, the first-order seismic wave equation is:

[0033]

[0034]

[0035]

[0036] Among them: S (x, y, z) is the source wave field, x, y, z are spatial coordinates, represents phase, t represents time, v(x,y,z) represents velocity model, represents the seismic wavelet at the position x0, y0, z0, p x (x,y,z),p y (x,y,z),p z (x, y, z) represent the stress field variables in the three spatial directions of x, y, and z respectively;

[0037] The second-order seismic wave field equation is:

[0038]

[0039] Preferably, the specific process of completing the forward continuation of the source wavefield of the optimized time for any shot seismic data in the shot gather data of step S3 is as follows:

[0040] In step S2, the seismic wave equation is added at x0, y0, z0 to generate the source wave field. Within the imaging space of the shot data, the seismic wave equation calculation formula in step S2 is used to optimize the forward extension of the source wave field along the direction of increasing time, that is, starting from time 0 and propagating to time T, where T is the propagation time length of the optimized source wave field. The source wave field at time 0-T is first compressed and then stored on the local disk.

[0041] Preferably, the specific process of completing the reverse time extension of the wave field of the detection point with optimized time for any shot seismic data in the S3 shot gather data is as follows:

[0042] In this step, the actual detection point position x, y, z1(x, y) of the seismic wave equation is added to the shot gather data of step S1. Within the imaging space of the shot gather data, the seismic wave equation calculation formula in this step is used to perform reverse time extension of the detection point wave field of the total trace length of the shot gather data along the direction from large to small recording time of the shot gather data of step S1, that is, starting from the maximum trace length time and continuing to reverse time extension to time 0;

[0043] ① For the first-order seismic wave field equation:

[0044]

[0045]

[0046]

[0047] Where: R(x,y,z) is the wave field of the receiver point; x,y,z are spatial coordinates; t represents time, the maximum time of which is the total length of the shot gather record in step 1; v(x,y,z) represents the velocity model; g(x,y,z1(x,y),t) represents the seismic record received by the receiver at the x,y,z1(x,y) position, where z1(x,y) represents the depth of the receiver that changes with the spatial (x,y) coordinates, including the changes in the position of the receiver on the same horizontal plane and on the undulating surface; px(x,y,z), p y (x,y,z),p z (x, y, z) represent the stress field variables in the three spatial directions of x, y, and z respectively;

[0048] ② For the second-order seismic wave field equation:

[0049]

[0050] Preferably, the propagation time length T of the source wavefield after the optimization of the source wavefield can be selected to be half or more of the shot gather data trace length in step 1, and less than the shot gather data trace length;

[0051] More precisely, it can be obtained according to the ray tracing method or the travel time calculation method based on the eikonal equation, specifically according to the maximum travel time value from the shot point to the imaging aperture range and depth range.

[0052] Preferably, in step S4, at any time step of the optimized time in the previous step, an imaging condition is applied to any imaging point position in the imaging space of the shot data, and the seismic wave field after the imaging condition is applied is accumulated within the 0-T time. The specific process of completing the reverse time migration imaging processing of the shot data in the shot gather data is as follows:

[0053] 31) The wave field of the detection point is reverse-time extended (in the direction of decreasing time) starting from the time when the total trace length of the shot gather data recorded in step S1 is the maximum. When the time Tt of the reverse-time extended wave field of the detection point is greater than T, no processing is performed. When the time Tt of the extended wave field of the detection point is equal to T, the source wave field corresponding to time T is first decompressed from the local disk, and the reverse-time imaging condition is applied to obtain the reverse-time imaging result of the shot data at time T.

[0054] 32) Then, when the wave field extension time of the detection point reaches time Tt, the source wave field at time Tt is decompressed from the local disk, and the reverse time imaging condition is applied to obtain the reverse time imaging result of the shot data at time Tt, thereby completing the reverse time imaging results of all times 0-T; and at the same time step, the absolute value operation of the source wave field value is performed to obtain the amplitude energy result at that time;

[0055] 34) Accumulate the reverse time imaging results at all 0-T moments in step 32) to obtain the reverse time migration result of the shot data; accumulate the amplitude energy results at all 0-T moments in step 32) to obtain the amplitude energy data volume of the shot data;

[0056] Preferably, the reverse time imaging condition applied at any time step is that the source wavefield and the receiver wavefield are calculated at the same time. Specific reverse time imaging conditions include:

[0057] a. Cross-correlation imaging conditions make the imaging results more stable;

[0058] b. Normalizing the reverse time imaging conditions to make the spatial amplitude energy distribution of the imaging results more uniform, with more prominent formation details and higher resolution;

[0059] c. The reverse time imaging condition based on traveling wave separation suppresses low-wavenumber noise interference during the imaging process, making the signal-to-noise ratio of the imaging results higher.

[0060] Preferably, the reverse time migration data volume in step S6 is optimized, and the optimization method includes:

[0061] 41) Perform amplitude compensation processing;

[0062] 42) Perform low wave number noise attenuation processing;

[0063] 43) According to the application requirements of the wave group characteristics, the processing result can also be multiplied by -1 to achieve a 180-degree phase rotation, which is used to eliminate the influence of the change in the phase of the input source wavelet in step 4 and the inability to consider the polarity of the output seismic wavelet after imaging, and obtain the final amplitude-compensated reverse-time imaging data volume.

[0064] Preferably, the amplitude compensation processing method is obtained by dividing the reverse time migration data volume of the entire working area in step S6 by the amplitude energy data volume;

[0065] The amplitude energy data body used as the denominator in the amplitude compensation process needs to be counted for its absolute value to obtain the maximum absolute value of the data body, and then a threshold value is set to limit the value range of the amplitude energy data body to obtain a modified amplitude energy data body;

[0066] Preferably, the method of setting a threshold value to limit the numerical range of the amplitude energy data volume includes:

[0067] A value between [80% and 98%] of the maximum absolute value of the amplitude energy data body is used as the maximum value selected for the modified amplitude energy data body, that is, a value in the amplitude energy data body that is greater than the maximum value selected is set as the maximum value;

[0068] A value between [1% and 10%] of the maximum absolute value of the amplitude energy data body is used as the minimum value of the modified amplitude energy data body, that is, a value in the amplitude energy data body that is less than the minimum value is set to the minimum value;

[0069] The above processing can eliminate the abnormal processing results caused by the denominator being zero or the value being too large or too small when dividing the global reverse time migration data volume by the global amplitude energy data volume.

[0070] The present invention also provides a phase-controllable and efficient reverse time migration imaging device, comprising:

[0071] The first acquisition unit is used to obtain the shot gather data and high-precision depth domain velocity model after fidelity processing of the seismic work area;

[0072] an adjustment unit for adjusting the phase of an input seismic wavelet for generating a source wavefield according to the type of seismic wave equation adopted;

[0073] The execution unit completes the forward extension of the source wave field with optimized time and the reverse extension of the receiver wave field with long total trace time for any shot seismic data in the shot gather data;

[0074] an accumulation unit, configured to apply an imaging condition to any imaging point in the imaging space of the shot data at any time step of the optimized time in step S3, and accumulate the seismic wavefield after the imaging condition is applied within the 0-T time to complete the reverse time migration imaging processing of the shot data in the shot gather data;

[0075] The second acquisition unit obtains a single-shot reverse-time imaging data volume and an amplitude energy data volume based on reverse-time migration imaging processing of any single-shot data in the completed shot gather data;

[0076] The first processing unit performs reverse time migration imaging processing on any single shot data in the shot gather data to obtain a reverse time imaging data volume and an amplitude energy data volume of the single shot;

[0077] The second processing unit is used to complete the reverse time migration processing of all shot data in the work area in step S5, and obtain the reverse time migration data volume and amplitude energy data volume of the entire work area by cutting and superimposing all single-shot reverse time imaging data volumes and amplitude energy data volumes, and complete the reverse time migration optimization processing.

[0078] Compared with the above-mentioned background technology, the present invention has the following beneficial effects: the phase-controllable and efficient reverse time migration imaging method of the present invention can solve the wave group difference between the reverse time migration imaging method and other existing imaging methods during the imaging process according to the type of seismic wave equation and the way of changing the source wavelet phase, and can also obtain reverse time migration results with approximately zero phase or arbitrary phase, effectively solving the wavelet space variation and wave field complexity problems existing in the post-stack phase shaping processing. At the same time, the reverse time imaging processing is more efficient and the processing results have the advantage of better spatial energy consistency, which is conducive to the accurate and detailed characterization of complex reservoirs or complex structures. BRIEF DESCRIPTION OF THE DRAWINGS

[0079] Figure 1 The initial depth domain velocity model of the embodiment of the present invention;

[0080] Figure 2 The depth domain velocity model after three rounds of fine tomographic inversion in an embodiment of the present invention;

[0081] Figure 3 The final depth domain imaging gathers and residual delay spectrum of the embodiment of the present invention;

[0082] Figure 4The final depth domain velocity model and imaging profile superimposed display diagram of the embodiment of the present invention;

[0083] Figure 5 A partially enlarged display of the zero-phase source seismic wavelet time-reverse imaging result according to an embodiment of the present invention;

[0084] Figure 6 A partially enlarged display of the Kirchhoff integration method prestack depth migration results according to an embodiment of the present invention;

[0085] Figure 7 A partially enlarged display of the reverse time imaging result of the seismic wavelet of the 90-degree phase source according to an embodiment of the present invention;

[0086] Figure 8 The result of reverse time imaging of seismic wavelet with 90-degree phase source according to the embodiment of the present invention;

[0087] Figure 9 Kirchhoff integration method prestack depth migration result diagram of the study area of ​​the embodiment of the present invention;

[0088] Figure 10 This is a flowchart of the phase-controllable and efficient reverse time migration imaging method of the present invention. DETAILED DESCRIPTION

[0089] The present invention will be further described below in conjunction with the embodiments and accompanying drawings:

[0090] like Figure 10 As shown, the phase-controllable and efficient reverse time migration imaging method includes the following steps:

[0091] 1. Obtain the high-fidelity shot gather data and high-precision depth-domain velocity model for the seismic area to be processed.

[0092] Shot gather data undergoes surface-consistent deconvolution and surface-consistent amplitude-energy compensation to eliminate the combined effects of source excitation (eliminating energy unevenness between shots) and receiver reception (eliminating wavefield differences between receivers), ensuring spatial consistency of seismic wavelet amplitude energy, phase, and frequency across all shot and receiver points. Both surface-consistent deconvolution and surface-consistent energy compensation utilize a long time window, encompassing at least 1 second of the target interval, to ensure that the processed seismic data retains the complete wavefield propagation effects and energy relationships, including seismic reflection and transmission.

[0093] For the velocity model, it is necessary to adopt a sophisticated three-dimensional velocity inversion method, including the full waveform inversion method, the multi-information constrained grid tomography velocity inversion method, etc.

[0094] The multi-information constrained grid tomographic velocity inversion method includes:

[0095] 11) Use prestack time migration velocity to perform time-to-depth conversion, and use VSP velocity or logging velocity, interpreted horizon information, and other information to constrain the construction of an initial depth domain velocity model. Complete the first round of prestack depth migration, output the first round of depth domain imaging volume and depth domain imaging gathers, and use a migration grid of 10 bins or less for the main and tie lines.

[0096] 12) Complete the residual delay spectrum calculation and seed point picking of the first round of depth domain imaging gathers, with the picking density in the depth direction being less than or equal to 40 times the depth step; and complete the three-dimensional grid tomographic velocity inversion, with the spatial direction inversion grid being less than or equal to 20 bins, and the depth direction inversion grid being less than or equal to 40 times the depth step. The three-dimensional grid tomographic inversion velocity model in this step is added to the initial depth domain model in step 11) to obtain a velocity model representing large-scale structural features.

[0097] 13) Complete the second round of pre-stack depth migration and output the second round of depth domain imaging volume and depth domain imaging gathers. The migration grid of the main survey line and the tie survey line should be less than or equal to 5 bins.

[0098] 14) Completing the residual delay spectrum calculation and seed point picking for the second round of depth domain imaging gathers, with the picking density in the depth direction being less than or equal to 20 times the depth step; and completing the three-dimensional grid tomographic velocity inversion, with the spatial inversion grid being less than or equal to 10 bins and the depth inversion grid being less than or equal to 20 times the depth step. The three-dimensional grid tomographic inversion velocity model in this step is added to the velocity model representing large-scale structural features in step 12), thereby obtaining a velocity model representing medium- and large-scale complex structures;

[0099] 15) Complete the third round of pre-stack depth migration and output the third round of depth domain imaging volume and depth domain imaging gathers. The migration grid of the main survey line and the tie survey line is less than or equal to 2 bins.

[0100] 16) Completing the residual delay spectrum calculation and seed point picking for the third round of depth domain imaging gathers, with the picking density in the depth direction being less than or equal to 20 times the depth step; and completing the three-dimensional grid tomographic velocity inversion, with the spatial inversion grid being less than or equal to 10 bins and the depth inversion grid being less than or equal to 20 times the depth step. The three-dimensional grid tomographic inversion velocity model in this step is added to the velocity model representing the large-scale complex structures in step 14) to obtain a velocity model representing the fine structure;

[0101] The number of iterations of the depth-domain velocity model is determined by whether the corresponding depth-domain imaging data volume meets the geological requirements (accurate structural imaging position, clear geological body imaging, clear fault and fracture characterization, etc.), the flattening degree of the imaging gather (the gather is flattened), and whether the energy cluster on the residual delay spectrum is zeroed. If the above standards are not met, the number of iterations of the depth-domain velocity model is increased to further improve the accuracy of the velocity model. Until a velocity model that characterizes fine structures that meets the standards is obtained, the iteration of the depth-domain velocity model ends.

[0102] The conditions for the structure matching between the velocity model and the imaging results are: the top and bottom depth positions of the velocity changes in the depth domain velocity model are consistent with the depth position and trend of the marker layer in the seismic imaging results; the local velocity changes in the depth domain velocity model are consistent with the local geological body change range in the seismic imaging results.

[0103] 2. According to the type of seismic wave equation used, adjust the phase of the input seismic wavelet that generates the source wavefield; the specific method includes the following steps:

[0104] The seismic wavelet for generating the source wavefield includes a zero-phase Ricker wavelet and a zero-phase frequency-domain band-limited seismic wavelet.

[0105] 1) To make the wave group characteristics of the reverse time imaging results consistent with those of the Kirchhoff depth imaging results, it is necessary to perform a 90-degree phase rotation on the input seismic wavelet that generates the source wavefield;

[0106] 2) If the reverse time imaging result is to have a wave group characteristic close to zero phase, the input seismic wavelet generating the source wavefield is not subjected to phase rotation;

[0107] 3) If you want the reverse time imaging result to have The wave group characteristics of the phase can be obtained by Phase rotation is achieved;

[0108] Preferably, the types of seismic wave equations used include first-order seismic wave equations and second-order seismic wave equations;

[0109] Taking the three-dimensional acoustic wave equation as an example, the first-order seismic wave equation is:

[0110]

[0111]

[0112]

[0113] Among them: S (x, y, z) is the source wave field, x, y, z are spatial coordinates, represents the phase, t represents the propagation time, v(x,y,z) represents the velocity model, represents the seismic wavelet at the position x0, y0, z0, p x (x,y,z),p y (x,y,z),p z (x, y, z) represent the stress field variables in the three spatial directions of x, y, and z respectively;

[0114] The second-order seismic wave field equation is:

[0115]

[0116] 3. For any shot seismic data in the shot gather data in step S1, complete the forward extension of the source wave field with optimized time and the reverse extension of the receiver wave field with long total trace time. The specific process is as follows:

[0117] a. The specific process of completing the forward extension of the source wave field with optimized time is as follows:

[0118] In step S2, the seismic wave equation is added at x0, y0, z0 to generate the source wave field. Within the imaging space of the shot data, the seismic wave equation calculation formula in step S2 is used to optimize the forward extension of the source wave field along the direction of increasing time, that is, starting from time 0 and propagating to time T, where T is the propagation time length of the source wave field after optimization. The source wave field at time 0-T is first compressed and then stored on the local disk.

[0119] b. The specific process of completing the reverse time extension of the detection point wave field with optimized time is as follows:

[0120] In this step, the actual detection point position x, y, z1(x, y) of the seismic wave equation is added to the shot gather data of step S1. Within the imaging space of the shot gather data, the seismic wave equation calculation formula in this step is used to perform reverse time extension of the detection point wave field of the total trace length of the shot gather data along the direction from large to small recording time of the shot gather data of step S1, that is, starting from the maximum trace length time and continuing to reverse time extension to time 0;

[0121] ① For the first-order seismic wave field equation:

[0122]

[0123]

[0124]

[0125] Where: R(x,y,z) is the wave field of the receiver point; x,y,z are spatial coordinates; t represents time, the maximum time of which is the total length of the shot gather record in step 1; v(x,y,z) represents the velocity model; g(x,y,z1(x,y),t) represents the seismic record received by the receiver at the x,y,z1(x,y) position, where z1(x,y) represents the depth of the receiver that changes with the spatial (x,y) coordinates, including the changes in the position of the receiver on the same horizontal plane and on the undulating surface; px(x,y,z), p y (x,y,z),p z (x,y,z) represents the stress field variables in the three spatial directions of x, y, and z respectively.

[0126] ② For the second-order seismic wave field equation:

[0127]

[0128] The length T of the source wavefield propagation time after the optimization processing can be selected to be half or more of the shot data track length in step 1, and less than the shot data track length; more precisely, it can be obtained according to the ray tracing method or the travel time calculation method based on the eikonal equation, specifically according to the maximum travel time value from the shot point to the imaging aperture range and depth range.

[0129] 4. At any time step optimized in step S3, the imaging condition is applied to any imaging point in the imaging space of the shot data, and the seismic wave field after the imaging condition is applied is accumulated within the 0-T time to complete the reverse time migration imaging processing of the shot data in the shot gather data. The specific process is as follows:

[0130] 1) The wave field of the detection point is reverse-time extended (towards a decreasing time direction) starting from the time when the total trace length of the shot gather data recorded in step S1 is the maximum. When the time Tt of the reverse-time extended wave field of the detection point is greater than T, no processing is performed. When the time Tt of the reverse-time extended wave field of the detection point is equal to T, the source wave field corresponding to time T is first decompressed from the local disk, and the reverse-time imaging condition is applied to obtain the reverse-time imaging result of the shot data at time T.

[0131] 2) Then, when the wavefield extension of the receiver point reaches time Tt, the source wavefield at time Tt is decompressed from the local disk and the reverse imaging condition is applied to obtain the reverse imaging result of the shot data at time Tt, thereby completing the reverse imaging results of all times 0-T. At the same time step, the absolute value of the source wavefield value is calculated to obtain the amplitude energy result at that time.

[0132] 3) Accumulating the reverse time imaging results at all 0-T moments in step 2) to obtain the reverse time migration results of the shot data; accumulating the amplitude energy results at all 0-T moments in step 2) to obtain the amplitude energy data volume of the shot data; the reverse time imaging conditions applied at any time step, i.e., the source wavefield and the receiver wavefield are calculated at the same time, specifically the reverse time imaging conditions include:

[0133] a. Cross-correlation imaging conditions make the imaging results more stable;

[0134] b. Normalizing the reverse time imaging conditions to make the spatial amplitude energy distribution of the imaging results more uniform, with more prominent formation details and higher resolution;

[0135] c. The reverse time imaging condition based on traveling wave separation suppresses low-wavenumber noise interference during the imaging process, making the signal-to-noise ratio of the imaging results higher.

[0136] 5. According to steps S3 and S4, the reverse time migration imaging processing of any single shot data in the shot gather data in step S1 is completed to obtain the single shot reverse time imaging data volume and amplitude energy data volume.

[0137] 6. According to step S5, the reverse time migration processing of all shot data in the work area is completed. By cutting and superimposing all single-shot reverse time imaging data volumes and amplitude energy data volumes, the reverse time migration data volume and amplitude energy data volume of the entire work area are obtained, and the reverse time migration optimization processing is completed.

[0138] The method for optimizing the reverse time migration data volume includes:

[0139] 1) Perform amplitude compensation processing;

[0140] 2) Perform low wave number noise attenuation processing;

[0141] 3) According to the application requirements of the wave group characteristics, the processing result can also be multiplied by -1 to achieve a 180-degree phase rotation, which is used to eliminate the influence of the phase change of the input source wavelet in step 4 and the inability to consider the polarity of the output seismic wavelet after imaging, and obtain the final amplitude-compensated reverse-time imaging data volume.

[0142] The amplitude compensation processing method is to obtain the amplitude energy data volume by dividing the reverse time migration data volume of the entire work area in step S6; the amplitude energy data volume used as the denominator in the amplitude compensation processing needs to be statistically analyzed for the absolute value of the amplitude energy data volume to obtain the maximum absolute value of the data volume, and then a threshold value is set to limit the value range of the amplitude energy data volume to obtain a modified amplitude energy data volume;

[0143] The method of setting a threshold value to limit the numerical range of the amplitude energy data body includes:

[0144] A value between [80% and 98%] of the maximum absolute value of the amplitude energy data body is used as the maximum value selected for the modified amplitude energy data body, that is, a value in the amplitude energy data body that is greater than the maximum value selected is set as the maximum value;

[0145] A value between [1% and 10%] of the maximum absolute value of the amplitude energy data body is used as the minimum value of the modified amplitude energy data body, that is, a value in the amplitude energy data body that is less than the minimum value is set to the minimum value;

[0146] The above processing can eliminate the abnormal processing results caused by the denominator being zero or the value being too large or too small when dividing the global reverse time migration data volume by the global amplitude energy data volume.

[0147] The present invention provides a phase-controllable and efficient reverse time migration imaging device, comprising:

[0148] The first acquisition unit is used to obtain the shot gather data and high-precision depth domain velocity model after fidelity processing of the seismic work area;

[0149] an adjustment unit for adjusting the phase of an input seismic wavelet for generating a source wavefield according to the type of seismic wave equation adopted;

[0150] The execution unit completes the forward extension of the source wave field with optimized time and the reverse extension of the receiver wave field with long total trace time for any shot seismic data in the shot gather data;

[0151] an accumulation unit, configured to apply an imaging condition to any imaging point in the imaging space of the shot data at any time step of the optimized time in step S3, and accumulate the seismic wavefield after the imaging condition is applied within the 0-T time to complete the reverse time migration imaging processing of the shot data in the shot gather data;

[0152] The second acquisition unit obtains a single-shot reverse-time imaging data volume and an amplitude energy data volume based on reverse-time migration imaging processing of any single-shot data in the completed shot gather data;

[0153] The first processing unit performs reverse time migration imaging processing on any single shot data in the shot gather data to obtain a reverse time imaging data volume and an amplitude energy data volume of the single shot;

[0154] The second processing unit is used to complete the reverse time migration processing of all shot data in the work area in step S5, and obtain the reverse time migration data volume and amplitude energy data volume of the entire work area by cutting and superimposing all single-shot reverse time imaging data volumes and amplitude energy data volumes, and complete the reverse time migration optimization processing.

[0155] Example 1

[0156] To make the purpose, technical solutions and advantages of the present invention clearer, the following will take the time-reversal imaging of carbonate reservoirs in a certain work area of ​​the Tadong exploration area of ​​Daqing Oilfield as an example and further describe the embodiments of the present invention in detail with reference to the accompanying drawings.

[0157] As attached Figure 10 As shown, the phase-controllable and efficient reverse time migration imaging method includes the following steps:

[0158] 1. Obtain the pre-processed shot gather data and high-precision depth-domain velocity model for the seismic work area to be processed. The bin size of the work area is 20 meters × 20 meters, and the migration depth step is 5 meters.

[0159] The shot gather data after fidelity preprocessing have undergone surface-consistent deconvolution and surface-consistent amplitude-energy compensation. The analysis window and application window include the target layer segment for 3 seconds (the total seismic recording trace time is 6 seconds). This makes the seismic wavelet amplitude energy, phase, and frequency of all shot points and receiver points spatially consistent. The processed seismic data retains the complete wave field propagation effects and energy relationships such as seismic reflection and transmission.

[0160] The high-precision velocity model is obtained by using a multi-information constrained grid tomographic velocity inversion method, which uses three rounds of velocity iteration from coarse to fine. The specific steps are as follows:

[0161] 1) Using pre-stack time migration velocity to perform time-depth conversion to construct the initial depth domain velocity model (e.g. Figure 1 The first round of prestack depth migration is completed, and the first round of depth domain imaging volume and depth domain imaging gathers are output. The migration grid of the main survey line and the tie survey line is 10 bins.

[0162] 2) Complete the residual delay spectrum calculation and seed point picking for the first round of depth-domain imaging gathers, with a depth-wise picking density of 200 meters (i.e., 40 times the depth step size); and complete the 3D grid tomographic velocity inversion, with a spatial inversion grid of 400 meters (i.e., 20 times the bin size) and a depth-wise inversion grid of 200 meters (i.e., 40 times the depth step size). The 3D grid tomographic inversion velocity model in this step is added to the initial depth-domain model in step 1) to obtain a velocity model representing large-scale structural features.

[0163] 3) Complete the second round of pre-stack depth migration and output the second round of depth domain imaging volume and depth domain imaging gathers. The migration grid for the main survey line and the tie survey line is 5 bins.

[0164] 4) Complete the residual delay spectrum calculation and seed point picking of the second round of depth domain imaging gathers, with a picking density of 100 meters in the depth direction (i.e., 20 times the depth step); and complete the three-dimensional grid tomographic velocity inversion, with a spatial inversion grid of 200 meters (i.e., 10 times the bin size) and a depth inversion grid of 100 meters (i.e., 20 times the depth step). The three-dimensional grid tomographic inversion velocity model in this step is added to the velocity model representing large-scale structural features in step 2) to obtain a velocity model representing complex structures at medium and large scales.

[0165] 5) Complete the third round of pre-stack depth migration and output the third round of depth domain imaging volume and depth domain imaging gathers. The migration grid of the main survey line and the tie survey line is less than or equal to 2 bins.

[0166] 6) Complete the calculation of the residual delay spectrum and seed point picking of the third round of depth domain imaging gathers. Considering the low signal-to-noise ratio of the data in the study area, the picking density in the depth direction is 100 meters (i.e., 20 times the depth step); and complete the three-dimensional grid tomographic velocity inversion. The spatial direction inversion grid is 200 meters (i.e., 10 times the bin size), and the depth direction inversion grid is 100 meters (i.e., 20 times the depth step). The three-dimensional grid tomographic inversion velocity model in this step is added to the velocity model representing the large-scale complex structure in step 4) to obtain the velocity model representing the fine structure (such as Figure 2 At this time, the depth domain imaging gather is flattened, and the energy cluster on its residual delay spectrum returns to zero (as shown in Figure 3 As shown in the figure), the complex structure imaging position is accurate, the geological body imaging is clear, the faults and fractures are clearly depicted, and the final depth domain velocity model is consistent with the target line depth domain imaging result structural trend (as shown in the figure). Figure 4 shown).

[0167] 2. In order to study the influence of the source wavefield seismic wavelet on the wave group phase characteristics of the reverse-time migration results and its difference from the conventional Kirchhoff depth domain seismic imaging method, the second-order seismic wave equation type was selected, and zero-phase seismic wavelet and 90-degree phase-rotated seismic wavelet were used as the source seismic wavelet to carry out reverse-time imaging processing application. The use of zero-phase seismic wavelet can make the reverse-time imaging results have a wave group characteristic of approximately zero phase, while the use of 90-degree phase-rotated seismic wavelet can make the reverse-time imaging results consistent with the wave group characteristics of the Kirchhoff depth imaging results.

[0168] The three-dimensional acoustic wave equation in the second-order form is as follows:

[0169]

[0170] Among them: S (x, y, z) is the source wave field, x, y, z are spatial coordinates, represents the phase, t represents the propagation time, v(x,y,z) represents the velocity model, represents the seismic wavelet at the position x0, y0, z0;

[0171] 3. In step 2, add the seismic wavelet that generates the source wave field at x0, y0, z0 in the seismic wave equation Within the imaging space of the shot data, the seismic wave equation calculation formula in step 2 is used to optimize the forward extension of the source wave field along the direction of increasing time. That is, starting from time 0, it propagates to T = 4s, which is 2 / 3 of the total trace time of 6s. The source wave field from 0 to 4s is first compressed and then stored on the local disk.

[0172] 4. Use the following second-order seismic wave equation to perform reverse time extension on any shot data in the shot set in step 1:

[0173]

[0174] Where: R(x,y,z) is the wave field of the receiver point; x,y,z are spatial coordinates; t represents time, with the maximum total trace length of the shot gather record in time step 1; v(x,y,z) represents the velocity model; g(x,y,z1(x,y),t) represents the seismic record received by the receiver at the x,y,z1(x,y) position, where z1(x,y) represents the depth of the receiver that changes with the spatial (x,y) coordinates, including changes in the receiver position on the same horizontal plane and on undulating surfaces.

[0175] That is, each sampling point of each seismic record in the shot set is added to the actual detection point position x, y, z1(x, y) according to the corresponding reverse time propagation time. Within the imaging space of the shot data, the seismic wave equation calculation formula in this step is used to reversely extend the detection point wave field along the direction from large to small total time of this step, that is, starting from 6s and continuing to reverse time extension to 0s.

[0176] 5. When the wave field of the detection point is reverse-time extended from 6s to 0s in step 3, if the reverse-time extended time Tt of the wave field of the detection point is greater than 4s, no processing is performed; when the extended time Tt of the wave field of the detection point is equal to 4s, the source wave field corresponding to the time of 4s is decompressed from the local disk and obtained, and the imaging value after the reverse-time imaging condition of the correlation method is multiplied by the source wave field and the wave field of the detection point is obtained, and the reverse-time imaging result of the shot data at the time of 4s is obtained; and so on. When the extended time of the wave field of the detection point reaches the time Tt, the source wave field at the time of Tt is decompressed from the local disk. The source wave field is obtained by multiplying the source wave field and the receiver wave field to obtain the imaging value after the reverse time imaging condition of the correlation method is applied, and the reverse time imaging result of the shot data at time Tt is obtained, thereby completing the reverse time imaging results of all times 0-4s; and at the same time step, the absolute value operation of the source wave field value is performed to obtain the amplitude energy result at time Tt; according to the above steps, the reverse time imaging results of all times 0-4s are completed and accumulated to obtain the reverse time migration result of the shot data; the amplitude energy results of all times 0-4s are completed and accumulated to obtain the amplitude energy data volume of the shot data;

[0177] 6. Following steps 2-5, the reverse time migration (RTM) processing of the 26,000-shot data for the entire work area was completed. The single-shot amplitude energy data volume and the single-shot RTM data volume were excised using the same excision function. The all-shot amplitude energy data volume and the all-shot RTM data volume were then superimposed, resulting in the full-area RTM data volume and the full-area amplitude energy data volume.

[0178] 7. Optimize the full-area RTM data volume in step 6. First, use amplitude compensation to divide the full-area RTM data volume in step 6 by the full-area amplitude energy data volume. The specific process is as follows:

[0179] The amplitude compensation process needs to first count the absolute value range of the amplitude energy data body of the entire area, thereby obtaining the maximum absolute value of the data body, and the minimum value is usually 0, and the value range of the amplitude energy data body of the entire area is adjusted to [5%, 90%] times the maximum absolute value, thereby obtaining the modified amplitude energy data body. On this basis, the diffusion filtering method is applied to perform low-wavenumber noise attenuation processing, and finally the zero-phase and 90-degree phase reverse-time imaging data bodies are obtained, whose imaging energy is more uniform and the signal-to-noise ratio is higher. From the analysis of the local reverse-time migration stacked profile, it can be seen that the reverse-time imaging results of the 90-degree phase-adjusted source seismic wavelet are consistent with the wave group phase characteristics of the Kirchhoff integral method imaging results, while the reverse-time imaging results of the zero-phase seismic wavelet have the characteristics of approximately zero phase (comparative analysis Figure 5 、 Figure 6 、 Figure 7 Finally, the reverse time imaging results of the work area using the 90-degree phase adjustment of the source seismic wavelet were selected ( Figure 8) and Kirchhoff integral method imaging results ( Figure 9 ) are consistent with the structural trends and wave group phase characteristics of the study area. The reverse time migration results can more clearly depict the stratigraphic details, have a higher signal-to-noise ratio, and more uniform transverse and vertical energy, which is helpful for reservoir prediction in the study area under the same seismic wave group characteristic conditions.

Claims

1. A phase-controllable and efficient reverse time migration imaging method, characterized by: The following steps are involved: S1. Obtain shot gather data and high-precision depth-domain velocity models after fidelity processing in the seismic area; S2. adjusting the phase of the input seismic wavelet that generates the source wavefield according to the type of seismic wave equation used; S3. For any shot seismic data in the shot data in step S1, the forward extension of the source wave field of the optimized time and the reverse extension of the total track long-time detection point wave field are completed; S4. At any time step optimized in step S3, applying the imaging condition to any imaging point in the imaging space of the shot gather data, and accumulating the seismic wavefield after the imaging condition is applied within the 0-T time, thereby completing the reverse time migration imaging process of the shot gather data in the shot gather data; S5. Complete the reverse time migration imaging processing of any single shot data in the shot gather data in step S1 according to step S3 and step S4 to obtain a single shot reverse time imaging data volume and an amplitude energy data volume; S6. Complete the reverse time migration processing of all shot data in the work area according to step S5, and obtain the reverse time migration data volume and amplitude energy data volume of the entire work area by cutting and superimposing all single-shot reverse time imaging data volumes and amplitude energy data volumes, and complete the reverse time migration optimization processing.

2. The phase-controllable and efficient reverse time migration imaging method according to claim 1, characterized in that: The shot gather data in step S1 needs to be subjected to surface-consistent deconvolution processing and surface-consistent amplitude energy compensation processing to ensure that the seismic wavelet amplitude energy, phase, and frequency of all shot points and receiver points are spatially consistent; The surface consistent deconvolution processing and surface consistent energy compensation processing both use a relatively long time window, at least including the target layer segment for more than 1 second, to ensure that the processed seismic data can retain the complete wave field propagation effect and energy relationship of seismic reflection and transmission.

3. The phase-controllable and efficient reverse time migration imaging method according to claim 2, characterized in that: The method for obtaining the high-precision depth-domain velocity model adopts a sophisticated three-dimensional velocity inversion method, which specifically includes a full waveform inversion method and a multi-information constrained grid tomography velocity inversion method.

4. The phase-controllable and efficient reverse time migration imaging method according to claim 3, characterized in that: The multi-information constrained grid tomographic velocity inversion method obtains a velocity model representing fine structures by iterating the depth-domain velocity model, and specifically includes the following steps: 11) Use prestack time migration velocity to perform time-to-depth conversion, and use VSP velocity or logging velocity and interpreted horizon information to constrain the construction of an initial depth-domain velocity model. Complete the first round of prestack depth migration, and output the first round of depth-domain imaging volume and depth-domain imaging gathers. The migration grid for the main and tie lines should be less than or equal to 10 bins. 12) Complete the residual delay spectrum calculation and seed point picking for the first round of depth-domain imaging gathers, with the depth-direction picking density being less than or equal to 40 times the depth step. Complete the three-dimensional grid tomographic velocity inversion, with the spatial inversion grid less than or equal to 20 bins and the depth-direction inversion grid less than or equal to 40 times the depth step. The three-dimensional grid tomographic inversion velocity model in this step is added to the initial depth-domain model from step 11) to obtain a velocity model representing large-scale structural features. 13) Complete the second round of pre-stack depth migration and output the second round of depth domain imaging volume and depth domain imaging gathers. The migration grid of the main survey line and the tie survey line should be less than or equal to 5 bins. 14) Complete the residual delay spectrum calculation and seed point picking for the second round of depth-domain imaging gathers, with the depth-wise picking density being less than or equal to 20 times the depth step. Complete the 3D grid tomographic velocity inversion, with the spatial inversion grid less than or equal to 10 bins and the depth-wise inversion grid less than or equal to 20 times the depth step. The 3D grid tomographic inversion velocity model in this step is added to the velocity model representing large-scale structural features in step 12) to obtain a velocity model representing complex structures at medium and large scales. 15) Complete the third round of pre-stack depth migration and output the third round of depth domain imaging volume and depth domain imaging gathers. The migration grid of the main survey line and the tie survey line is less than or equal to 2 bins. 16) Complete the residual delay spectrum calculation and seed point picking for the third round of depth-domain imaging gathers, with a depth-wise picking density of less than or equal to 20 times the depth step. Complete the 3D grid tomographic velocity inversion, with a spatial inversion grid of less than or equal to 10 bins and a depth-wise inversion grid of less than or equal to 20 times the depth step. The 3D grid tomographic inversion velocity model in this step is added to the velocity model representing medium- and large-scale complex structures in step 14) to obtain a velocity model representing fine structures. 17) When the depth-domain imaging data volume corresponding to the depth-domain velocity model meets the geological requirements, the velocity model matches the imaging result structure, and the imaging gather flattening standard, a velocity model that characterizes the fine structure that meets the standard is obtained; and the iteration of the depth-domain velocity model ends.

5. The phase-controllable and efficient reverse time migration imaging method according to claim 4, characterized in that: The factors that meet the geological requirements in step 17) are accurate structural imaging positions, clear geological body imaging, and clear fault and fracture depictions; The flattening standard requirement of the imaging gather is: the energy cluster on the residual delay spectrum of the depth domain imaging gather returns to zero; The conditions for the structure matching between the velocity model and the imaging results are: the top and bottom depth positions of the velocity changes in the depth domain velocity model are consistent with the depth position and trend of the marker layer in the seismic imaging results; the local velocity changes in the depth domain velocity model are consistent with the local geological body change range in the seismic imaging results.

6. The phase-controllable and efficient reverse time migration imaging method according to claim 1, characterized in that: The seismic wavelet for generating the source wavefield in step S2 includes a zero-phase Ricker wavelet and a zero-phase frequency-domain band-limited seismic wavelet.

7. The phase-controllable and efficient reverse time migration imaging method according to claim 6, characterized in that: The method for adjusting the phase of the input seismic wavelet for generating the source wavefield using the seismic wave equation type adopted in step S2 comprises the following steps: 21) To make the wave group characteristics of the reverse time imaging results consistent with those of the Kirchhoff depth imaging results, it is necessary to perform a 90-degree phase rotation on the input seismic wavelet that generates the source wavefield; 22) If the reverse time imaging result is to have a wave group characteristic close to zero phase, the input seismic wavelet generating the source wavefield should not be phase rotated; 23) To make the reverse time imaging result have The wave group characteristics of the phase are obtained by Phase rotation is achieved.

8. The phase-controllable and efficient reverse time migration imaging method according to claim 7, characterized in that: The types of seismic wave equations used include first-order seismic wave equations and second-order seismic wave equations; The first-order form of the seismic wave equation is: ， ， ， ； in: is the source wave field, is the spatial coordinate, Represents the phase, Represents time, represents the velocity model, represent The seismic wavelet at the location, 、 、 Respectively represent Stress field variables in three spatial directions; The second-order seismic wave field equation is: 。 9. The phase-controllable and efficient reverse time migration imaging method according to claim 1 or 6, characterized in that: The specific process of completing the forward continuation of the source wavefield of the optimized time for any shot seismic data in the shot gather data of step S3 is as follows: In step S2, the earthquake wave equation The seismic wavelet that generates the source wave field is added to the position Within the imaging space of the shot gather data, the seismic wave equation calculation formula in step S2 is used to optimize the forward extension of the source wave field along the direction of increasing time, that is, starting from time 0 and propagating to time T, where T is the propagation time length of the optimized source wave field. The source wave field at time 0-T is first compressed and then stored on the local disk.

10. The phase-controllable and efficient reverse time migration imaging method according to claim 1 or 6, characterized in that: The specific process of completing the reverse time extension of the wave field of the optimized detection point for any shot seismic data in the S3 shot gather data is as follows: The actual detection point position of the seismic wave equation in this step Add the shot gather data from step S1, and within the imaging space of the shot gather data, use the seismic wave equation calculation formula in this step to perform reverse time extension of the detection point wave field of the total trace length of the shot gather data along the direction of the recording time of the shot gather data from large to small, that is, starting from the maximum trace length time and continuing to reverse time extension to time 0; ① For the first-order seismic wave field equation: ， ， ， ； in: is the detection point wave field; is the spatial coordinate; represents the time, the maximum time of which is the total length of the shot gather recorded in step 1; represents the velocity model; represent The earthquake record received by the geophone at location The depth of the detector varies with the space The coordinates change, including the changes of the detection point position on the same horizontal plane and on the undulating surface; 、 、 Respectively represent Stress field variables in three spatial directions; ② For the second-order seismic wave field equation: 。 11. The phase-controllable and efficient reverse time migration imaging method according to claim 5, characterized in that: The optimized source wavefield propagation time length T is selected to be half or more of the shot gather data trace length in step 1 and less than the shot gather data trace length.

12. The phase-controllable and efficient reverse time migration imaging method according to claim 11, characterized in that: The propagation time length T of the earthquake source wave field is obtained by the ray tracing method or the travel time calculation method based on the eikonal equation.

13. The phase-controllable and efficient reverse time migration imaging method according to claim 11, characterized in that: The source wave field propagation time length T is obtained based on the maximum travel time value from the shot point to the imaging aperture range and depth range.

14. The phase-controllable and efficient reverse time migration imaging method according to claim 1, characterized in that: In the step S4, at any time step after the time is optimized in the previous step, the imaging condition is applied to any imaging point position in the imaging space of the shot gather data, and the seismic wave field after the imaging condition is applied is accumulated within the 0-T time. The specific process of completing the reverse time migration imaging processing of the shot gather data in the shot gather data is as follows: 31) The wave field of the detection point is reverse-time extended from the time corresponding to the maximum trace length of the shot gather data recorded in step S1. When the time Tt of the reverse-time extended wave field of the detection point is greater than T, no processing is performed. When the time Tt of the reverse-time extended wave field of the detection point is equal to T, the source wave field corresponding to time T is first decompressed from the local disk, and the reverse-time imaging condition is applied to obtain the reverse-time imaging result of the shot gather data at time T. 32) Then, when the wave field extension time of the receiver point reaches time Tt, the source wave field at time Tt is decompressed from the local disk, and the reverse time imaging condition is applied to obtain the reverse time imaging result of the shot gather data at time Tt, thereby completing the reverse time imaging results of all times 0-T; and at the same time step, the absolute value operation of the source wave field value is performed to obtain the amplitude energy result at that time; 34) Accumulate the reverse time imaging results of all 0-T moments in step 32) to obtain the reverse time migration result of the shot gather data; accumulate the amplitude energy results of all 0-T moments in step 32) to obtain the amplitude energy data volume of the shot gather data.

15. The phase-controllable and efficient reverse time migration imaging method according to claim 1, characterized in that: The reverse time migration data volume in step S6 is optimized. The optimization method includes: 1) Perform amplitude compensation processing; 2) Perform low wavenumber noise attenuation processing.

16. The phase-controllable and efficient reverse time migration imaging method according to claim 15, characterized in that: According to the application requirements of wave group characteristics, the low-wavenumber noise attenuation processing result is multiplied by -1 to achieve a 180-degree phase rotation, which is used to eliminate the influence of the phase change of the input source wavelet in step S4 and the inability to consider the polarity of the output seismic wavelet after imaging, and obtain the final amplitude-compensated reverse-time imaging data volume.

17. A phase-controllable and efficient reverse time migration imaging method according to claim 16, characterized in that The amplitude compensation processing method is obtained by dividing the reverse time migration data volume of the entire work area in step S6 with the amplitude energy data volume; The amplitude energy data body used as the denominator in the amplitude compensation process needs to be statistically analyzed for the absolute value of the amplitude energy data body to obtain the maximum absolute value of the data body, and then a threshold value is set to limit the value range of the amplitude energy data body to obtain a modified amplitude energy data body.

18. The phase-controllable and efficient reverse time migration imaging method according to claim 17, characterized in that: The method of setting a threshold value to limit the numerical range of the amplitude energy data body includes: A value between [80% and 98%] of the maximum absolute value of the amplitude energy data body is used as the maximum value selected for the modified amplitude energy data body, that is, a value in the amplitude energy data body that is greater than the maximum value selected is set as the maximum value; A value between [1% and 10%] of the maximum absolute value of the amplitude energy data body is used as the minimum value of the modified amplitude energy data body, that is, the value in the amplitude energy data body that is less than the minimum value is set to the minimum value; The above processing eliminates the abnormal processing results caused by the denominator being zero or the value being too large or too small when dividing the full-area reverse time migration data volume by the full-area amplitude energy data volume.

19. A phase-controllable and efficient reverse time migration imaging device, characterized by: The device realizes phase-controllable and efficient reverse time migration imaging by utilizing the method according to any one of claims 1 to 18, comprising: The first acquisition unit is used to obtain the shot gather data and high-precision depth domain velocity model after fidelity processing of the seismic work area; an adjustment unit for adjusting the phase of an input seismic wavelet for generating a source wavefield according to the type of seismic wave equation adopted; The execution unit completes the forward extension of the source wave field with optimized time and the reverse extension of the receiver wave field with long total trace time for any shot seismic data in the shot gather data; an accumulation unit, configured to apply an imaging condition to any imaging point in the imaging space of the shot gather data at any time step of the optimized time in step S3, and accumulate the seismic wave field after the imaging condition is applied within the 0-T time to complete the reverse time migration imaging processing of the shot gather data in the shot gather data; The second acquisition unit obtains a single-shot reverse-time imaging data volume and an amplitude energy data volume based on reverse-time migration imaging processing of any single-shot data in the completed shot gather data; The first processing unit performs reverse time migration imaging processing on any single shot data in the shot gather data to obtain a reverse time imaging data volume and an amplitude energy data volume of the single shot; The second processing unit is used to complete the reverse time migration processing of all shot data in the work area in step S5, and obtain the reverse time migration data volume and amplitude energy data volume of the entire work area by cutting and superimposing all single-shot reverse time imaging data volumes and amplitude energy data volumes, and complete the reverse time migration optimization processing.

Images