Phase-controllable efficient reverse time migration imaging method and device

Abstract

The invention provides a phase-controllable and efficient reverse time migration imaging method, and mainly solves the problems that the wavelet space variant and the complex wave field cannot be solved during post-stack phase shaping processing by the existing reverse time imaging method. The method comprises the following steps: S1, acquiring shot gather data and a high-precision depth domain velocity model after fidelity processing; s2, adjusting and generating the phase of an input seismic wavelet of a seismic source wave field; s3, for any shot of seismic data, completing optimization time seismic source wave field forward continuation and total trace length time detection point wave field reverse time continuation, applying imaging conditions, completing reverse time migration processing of any shot of seismic data, and obtaining a shot reverse time imaging data body and an amplitude energy data body; and S4, obtaining a reverse time migration data volume and an amplitude energy data volume of the whole work area. According to the method, the corresponding seismic source wavelet phase is changed according to the type of the adopted seismic wave equation, self-adaptive change of the imaging processing result phase is achieved in the reverse time migration process, and 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 technology is a technology that approximates seismic wave equations less, can solve multi-value travel time problems, is suitable for imaging problems in situations where any steep inclination angle and velocity changes dramatically in the horizontal and vertical directions, and can also accurately image multiple waves, rotating reflection waves, and other wave fields that are usually considered to be interference wave types. Compared with Kirchhoff prestack depth imaging technology and single-way wave seismic depth imaging technology, seismic wave time-reversal imaging technology has significant technical advantages and is one of the technologies with more mature theories and the most accurate imaging in the current high-precision seismic imaging technology series.

[0003] However, the problem with reverse time migration imaging is that the current reverse time migration imaging results have large differences in wave group phase characteristics from conventional Kirchhoff and other migration results. Some wave group characteristics are similar, and some 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 impossible to solve the problems of wave group characteristic differences introduced by sub-wave space variation and observation azimuth 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 a lot of 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 that the imaging stacking results in the background technology cannot solve the phase correction of the deep variation and space variation of the deep domain seismic wavelet, and provide a phase-controllable and efficient reverse time migration imaging method, which can change the corresponding source wavelet phase according to the type of seismic wave equation adopted, and efficiently realize the adaptive change of the phase of the imaging processing result in the reverse time migration process, thereby effectively solving the wavelet space variation and wave field complexity problems in the post-stack phase shaping process, 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 through 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 model after fidelity processing in the seismic work 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, complete the forward extension of the source wave field with optimized time and the reverse extension of the detection point wave field with a total time;

[0009] S4. At any time step of the optimized time in step S3, 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 to complete the reverse time migration imaging processing of the shot data in the shot gather data;

[0010] S5. According to step S3 and step 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;

[0011] S6. According to step S5, the reverse time migration processing of all shot data in the work area is completed, and the reverse time migration data volume and amplitude energy data volume of the whole work area are obtained by cutting and superimposing all single shot reverse time imaging data volumes and amplitude energy data volumes, and the reverse time migration optimization processing is completed.

[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 detection points are spatially consistent;

[0013] The surface consistent deconvolution processing and the 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 using 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 tomography 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 pre-stack time migration velocity to perform time-to-depth conversion, and use VSP velocity or logging velocity, interpretation layer information and other information to constrain and construct the initial depth domain velocity model, complete the first round of pre-stack depth migration, output the first round of depth domain imaging volume and depth domain imaging gathers, and the migration grid of the main survey line and the connecting survey line is less than or equal to 10 bins;

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

[0018] 13) Complete the second round of pre-stack depth migration, output the second round of depth domain imaging volume and depth domain imaging gathers, and the migration grid of the main survey line and the contact survey line is less than or equal to 5 bins;

[0019] 14) Complete the residual delay spectrum calculation and seed point picking of the second round of depth domain imaging gathers, and the picking density in the depth direction is less than or equal to 20 times the depth step; and complete the three-dimensional grid tomography velocity inversion, the spatial direction inversion grid is less than or equal to 10 facets, and the depth direction inversion grid is less than or equal to 20 times the depth step. The three-dimensional grid tomography inversion velocity model in this step is added to the velocity model characterizing the large-scale structural features in step 12), thereby obtaining a complex structural velocity model characterizing the medium and large scales;

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

[0021] 16) Complete the calculation of the residual delay spectrum and seed point picking of 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 complete the three-dimensional grid tomographic velocity inversion, with the inversion grid in the spatial direction being less than or equal to 10 facets, and the inversion grid in the depth direction 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 structure 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; 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 characterization;

[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 returns to zero.

[0025] The conditions for matching the velocity model with 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 type of seismic wave equation adopted in step S2, and the method for adjusting the phase of the input seismic wavelet for generating the source wave field, comprises the following steps:

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

[0029] 22) If the reverse time imaging result is to have a wave group characteristic of approximately zero phase, the input seismic wavelet that generates the source wave field is not phase rotated;

[0030] 23) 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;

[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 wave field 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 with the seismic wavelet generating the source wave field at x0, y0, z0. 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, and 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 data of step S1. Within the imaging space of the shot data, the calculation formula of the seismic wave equation in this step is used to perform reverse time extension of the detection point wave field of the total trace length of the shot data along the direction from large to small of the recording time of the shot data of step S1, that is, starting from the maximum trace length moment, reverse time extension has been performed to 0 moment;

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

[0044]

[0045]

[0046]

[0047] Where: R(x,y,z) is the wave field of the detection point; x,y,z are spatial coordinates; t represents time, and its maximum time 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 detector at the position x,y,z1(x,y), where z1(x,y) represents the depth of the detector that changes with the spatial (x,y) coordinates, including the changes in the detection point position 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 seismic source wave field after the seismic source wave field optimization 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 the depth range.

[0052] Preferably, in the step S4, at any time step of the optimized time in the previous step, the 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, and 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 reversely extended (towards the direction of decreasing time) from the maximum time of the total trace length of the shot data recorded in step S1. When the time Tt of the reverse extension of the wave field of the detection point is greater than T, no processing is performed; when the time Tt of the extension of the wave field of the detection point is equal to T, the source wave field corresponding to the time T is first decompressed from the local disk, and the reverse imaging condition is applied to obtain the reverse imaging result of the shot data at the time T;

[0054] 32) Then, when the time of the wave field extension of the detection point reaches the time Tt, the source wave field at the 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 the time Tt, thereby completing the reverse time imaging results of all the 0-T times; and at the same time step, the absolute value of the source wave field value is calculated to obtain the amplitude energy result at the time;

[0055] 34), cumulatively adding all the reverse time imaging results at the 0-T moment in step 32), to obtain the reverse time migration result of the shot data; cumulatively adding all the amplitude energy results at the 0-T moment in step 32), to obtain the amplitude energy data body of the shot data;

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

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

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

[0059] c. The reverse time imaging condition based on traveling wave separation suppresses the low-wavenumber noise interference during the imaging process, making the signal-to-noise ratio of the imaging result 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 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 of the phase of the input source wavelet in step 4 and the inability to consider the polarity of the output seismic wavelet after imaging, so as to obtain the reverse time imaging data body with the final amplitude compensation processing.

[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 with 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 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;

[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%, 98%] of the maximum absolute value of the amplitude energy data body is used as the maximum value selection value of the modified amplitude energy data body, that is, a value in the amplitude energy data body that is greater than the maximum value selection value is set as the maximum value;

[0068] A value between [1%, 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 result caused by the denominator being zero or the value being too large or too small when dividing the whole-area reverse time migration data body and the whole-area amplitude energy data body.

[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 after fidelity processing and the high-precision depth domain velocity model of the seismic work area;

[0072] An adjustment unit, for adjusting the phase of an input seismic wavelet for generating a source wave field 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 detection point wave field with long total trace time for any shot seismic data in the shot gather data;

[0074] An accumulation unit is used for applying imaging conditions at any imaging point position in the imaging space of the shot data at any time step of the optimization time in step S3, and accumulating the seismic wave field after the imaging conditions are applied within the 0-T time to complete the reverse time migration imaging processing of the shot data in the shot gather data;

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

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

[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 whole 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 approximate zero phase or arbitrary phase, effectively solving the wavelet space variation and wave field complexity problems existing in post-stack phase shaping processing. At the same time, the reverse time imaging processing is more efficient, and the processing result has 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 the 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 superposition 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 local enlarged display of the Kirchhoff integration method prestack depth migration result according to an embodiment of the present invention;

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

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

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

[0088] Fig.10 It is a flowchart of the phase-controllable and efficient reverse time migration imaging method of the present invention. Specific implementation method:

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

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

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

[0092] For the shot data, it is necessary to go through surface consistent deconvolution and surface consistent amplitude energy compensation to eliminate the combined influence of source excitation factors (eliminating energy unevenness between shots) and detector reception factors (eliminating wave field differences between detectors), and ensure that the amplitude energy, phase, and frequency of the seismic wavelets at all shot points and detector points are spatially consistent. Among them, surface consistent deconvolution and surface consistent energy compensation both use a longer time window, at least including the target layer segment for more than 1s, to ensure that the processed seismic data can retain the complete wave field propagation effects and energy relationships such as 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 tomography velocity inversion method includes:

[0095] 11) Use pre-stack time migration velocity to perform time-to-depth conversion, and use VSP velocity or logging velocity, interpretation layer information and other information to constrain and construct the initial depth domain velocity model, complete the first round of pre-stack depth migration, output the first round of depth domain imaging volume and depth domain imaging gathers, and the migration grid of the main survey line and the connecting survey line is less than or equal to 10 bins;

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

[0097] 13) Complete the second round of pre-stack depth migration, output the second round of depth domain imaging volume and depth domain imaging gathers, and the migration grid of the main survey line and the contact survey line is less than or equal to 5 bins;

[0098] 14) Complete the residual delay spectrum calculation and seed point picking of the second round of depth domain imaging gathers, and the picking density in the depth direction is less than or equal to 20 times the depth step; and complete the three-dimensional grid tomography velocity inversion, the spatial direction inversion grid is less than or equal to 10 facets, and the depth direction inversion grid is less than or equal to 20 times the depth step. The three-dimensional grid tomography inversion velocity model in this step is added to the velocity model characterizing the large-scale structural features in step 12), thereby obtaining a complex structural velocity model characterizing the medium and large scales;

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

[0100] 16) Complete the calculation of the residual delay spectrum and seed point picking of 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 complete the three-dimensional grid tomographic velocity inversion, with the inversion grid in the spatial direction being less than or equal to 10 facets, and the inversion grid in the depth direction 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 structure in step 14) to obtain a velocity model representing the fine structure;

[0101] Here, the number of iterations of the depth domain velocity model is determined based on 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 group 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 matching the velocity model with 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 adopted, the phase of the input seismic wavelet for generating the source wave field is adjusted; the specific method comprises 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) If the wave group characteristics of the reverse time imaging result are to be consistent with those of the Kirchhoff depth imaging result, it is necessary to perform a 90-degree phase rotation on the input seismic wavelet that generates the source wave field;

[0106] 2) If the reverse time imaging result is to have a wave group characteristic of approximately zero phase, the input seismic wavelet that generates the source wave field is not phase rotated;

[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 detection point wave field with 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 with the seismic wavelet generating the source wave field at x0, y0, z0. 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 the optimization of the source wave field, and 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 data of step S1. Within the imaging space of the shot data, the calculation formula of the seismic wave equation in this step is used to perform reverse time extension of the detection point wave field of the total trace length of the shot data along the direction from large to small of the recording time of the shot data of step S1, that is, starting from the maximum trace length moment, reverse time extension has been performed to 0 moment;

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

[0122]

[0123]

[0124]

[0125] Where: R(x,y,z) is the wave field of the detection point; x,y,z are spatial coordinates; t represents time, and its maximum time 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 detector at the position x,y,z1(x,y), where z1(x,y) represents the depth of the detector that changes with the spatial (x,y) coordinates, including the changes in the detection point position 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.

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

[0127]

[0128] The propagation time length T of the source wave field 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 the depth range.

[0129] 4. In step S3, at any time step of the optimized time, the 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 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 reversely extended (towards the direction of decreasing time) from the maximum total path length of the shot data recorded in step S1. When the time Tt of the reverse extension of the wave field of the detection point is greater than T, no processing is performed; when the time Tt of the extension of the wave field of the detection point is equal to T, the source wave field corresponding to the time T is first decompressed from the local disk, and the reverse imaging condition is applied to obtain the reverse imaging result of the shot data at the time T;

[0131] 2) Then, when the time of the wave field extension of the detection point reaches the time Tt, the source wave field at the 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 the time Tt, thereby completing the reverse time imaging results of all the 0-T times; and at the same time step, the absolute value of the source wave field value is calculated to obtain the amplitude energy result at that time;

[0132] 3) Accumulate the reverse time imaging results of all 0-T moments in step 2) to obtain the reverse time migration results of the shot data; Accumulate the amplitude energy results of 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, that is, the source wave field and the receiver wave field are calculated at the same time, and the specific reverse time imaging conditions include:

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

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

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

[0136] 5. According to step S3 and step 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, and the reverse time migration data volume and amplitude energy data volume of the whole work area are obtained by cutting and superimposing all single shot reverse time imaging data volumes and amplitude energy data volumes, 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 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 of the input source wavelet being changed in step 4 and the polarity of the output seismic wavelet after imaging cannot be considered, and obtain the reverse time imaging data body with the final amplitude compensation processing.

[0142] The amplitude compensation processing method is obtained by dividing the total work area reverse time migration data volume and the amplitude energy data volume in step S6; the amplitude energy data volume used as the denominator in the amplitude compensation processing needs to count the absolute value of the amplitude energy data volume to obtain the maximum absolute value of the data volume, and then set a threshold value to limit the value range of the amplitude energy data volume to obtain the 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%, 98%] of the maximum absolute value of the amplitude energy data body is used as the maximum value selection value of the modified amplitude energy data body, that is, a value in the amplitude energy data body that is greater than the maximum value selection value is set as the maximum value;

[0145] A value between [1%, 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 result caused by the denominator being zero or the value being too large or too small when dividing the whole-area reverse time migration data body and the whole-area amplitude energy data body.

[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 after fidelity processing and the high-precision depth domain velocity model of the seismic work area;

[0149] An adjustment unit, for adjusting the phase of an input seismic wavelet for generating a source wave field 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 detection point wave field with long total trace time for any shot seismic data in the shot gather data;

[0151] An accumulation unit is used for applying imaging conditions at any imaging point position in the imaging space of the shot data at any time step of the optimization time in step S3, and accumulating the seismic wave field after the imaging conditions are applied within the 0-T time to complete the reverse time migration imaging processing of the shot data in the shot gather data;

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

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

[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 whole 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] In order to make the purpose, technical scheme and advantages of the present invention clearer, the following will take the reverse time imaging of carbonate reservoirs in a certain work area of ​​the Tadong exploration area of ​​Daqing Oilfield as an example, and further describe the implementation mode of the present invention in detail in combination with the accompanying drawings.

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

[0158] 1. Obtain the shot gather data and high-precision depth domain velocity model after fidelity preprocessing in 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 data after fidelity preprocessing have been processed by surface consistent deconvolution and surface consistent amplitude energy compensation. The analysis window and application window include the target layer segment for 3s (the total trace time of seismic record is 6s), so that the amplitude energy, phase and frequency of seismic wavelets of all shot points and detection points are 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 adopts the multi-information constrained grid tomography velocity inversion method, which is obtained by three rounds of velocity iteration from coarse to fine. The specific steps are as follows:

[0161] 1) Use prestack time migration velocity to perform time-depth conversion to construct the initial depth domain velocity model (such as Figure 1 As shown in the figure, the first round of pre-stack 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 contact survey line is 10 facets.

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

[0163] 3) Complete the second round of pre-stack depth migration, output the second round of depth domain imaging volume and depth domain imaging gathers, and the migration grid of the main survey line and the contact 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 tomography velocity inversion, with a spatial direction inversion grid of 200 meters (i.e., 10 times the bin size), and a depth direction inversion grid of 100 meters (i.e., 20 times the depth step). The three-dimensional grid tomography inversion velocity model in this step is added to the velocity model characterizing large-scale structural features in step 2), thereby obtaining a complex structural velocity model characterizing medium and large scales;

[0165] 5) Complete the third round of pre-stack depth migration, output the third round of depth domain imaging volume and depth domain imaging gathers, and the migration grid of the main survey line and the contact 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 that the signal-to-noise ratio of the data in the study area is relatively low, the picking density in the depth direction is 100 meters (i.e., 20 times the depth step); and complete the three-dimensional grid tomography velocity inversion. The spatial direction inversion grid is 200 meters (i.e., 10 times the surface element size), and the depth direction inversion grid is 100 meters (i.e., 20 times the depth step). The three-dimensional grid tomography 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 At the same time, the complex structure imaging position is accurate, the geological body imaging is clear, the faults and fractures are clearly portrayed, and the final depth domain velocity model is consistent with the target line depth domain imaging result structural trend (as shown in Figure 1). 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 seismic wavelet after 90-degree phase rotation 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 result have a wave group characteristic of approximately zero phase, while the use of seismic wavelet after 90-degree phase rotation can make the reverse-time imaging result consistent with the wave group characteristics of the Kirchhoff depth imaging result.

[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 of 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, and the source wave field of 0-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 detection point; x,y,z are spatial coordinates; t represents time, and its maximum time step is the total length of the shot set 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 detector at the position x,y,z1(x,y), where z1(x,y) represents the depth of the detector that changes with the spatial (x,y) coordinates, including the changes in the detection point position on the same horizontal plane and on the undulating surface.

[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 in the total time of this step, that is, starting from 6s, and reversely extending to 0s.

[0176] 5. In the process of reverse time extension of the detection point wave field from 6s to 0s in step 3, when the reverse time extension time Tt of the detection point wave field is greater than 4s, no processing is performed; when the extension time Tt of the detection point wave field is equal to 4s, decompress and obtain the source wave field corresponding to the time of 4s from the local disk, and multiply the source wave field and the detection point wave field to obtain the imaging value after the reverse time imaging condition of the correlation method is applied, and obtain the reverse time imaging result of the shot data at the time of 4s; and so on, when the extension time of the detection point wave field reaches the time Tt, decompress and obtain the source wave field at the time Tt from the local disk The wave field is obtained by multiplying the source wave field and the detection point 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 all the reverse time imaging results at time 0-4s; and at the same time step, the absolute value operation is performed on the source wave field value to obtain the amplitude energy result at time Tt; according to the above steps, all the reverse time imaging results at time 0-4s are completed and accumulated to obtain the reverse time migration result of the shot data; all the amplitude energy results at time 0-4s are completed and accumulated to obtain the amplitude energy data body of the shot data;

[0177] 6. According to step 2 to step 5, the reverse time migration processing of 26,000 shot data in the whole work area is completed. The single shot amplitude energy data volume and the single shot reverse time migration data volume are cut by using the same cutting function, and then all shot amplitude energy data volumes and all shot reverse time migration data volumes are superimposed respectively, finally obtaining the reverse time migration data volume and the amplitude energy data volume of the whole area;

[0178] 7. The whole-area RTM data volume in step 6 is optimized by first using amplitude compensation processing to divide the whole-area RTM data volume in step 6 with the whole-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. And on this basis, the diffusion filtering method is applied to perform low-wavenumber noise attenuation processing, and finally the reverse-time imaging data body of zero phase and 90-degree phase is obtained, and its 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 result of the 90-degree phase-adjusted source seismic wavelet is consistent with the wave group phase characteristics of the Kirchhoff integral method imaging result, while the reverse-time imaging result of the zero-phase seismic wavelet has the characteristics of approximate zero phase (comparative analysis Figure 5 , Figure 6 , Figure 7 ). Finally, the reverse time imaging results of the work area with 90-degree phase adjustment of the source seismic wavelet were selected ( Figure 8) and Kirchhoff integral method imaging results ( Fig. 9 ) are consistent with the structural trends and wave group phase characteristics of the same seismic wave group. The reverse time migration results can depict the stratigraphic details more clearly, 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 in that: The following steps are involved: S1. Obtain shot gather data and high-precision depth domain velocity model after fidelity processing in the seismic work 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, complete the forward extension of the source wave field with optimized time and the reverse extension of the detection point wave field with a total time; S4. At any time step of the optimized time in step S3, 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 to complete the reverse time migration imaging processing of the shot data in the shot gather data; S5. According to step S3 and step 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; S6. According to step S5, the reverse time migration processing of all shot data in the work area is completed, and the reverse time migration data volume and amplitude energy data volume of the whole work area are obtained by cutting and superimposing all single shot reverse time imaging data volumes and amplitude energy data volumes, and the reverse time migration optimization processing is completed.

2. The phase-controllable and efficient reverse time migration imaging method according to claim 1, characterized in that: The shot data in step S1 needs to be processed by surface consistency deconvolution and surface consistency amplitude energy compensation 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 the 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 is obtained by using a sophisticated three-dimensional velocity inversion method, specifically including a full waveform inversion method and a multi-information constrained grid tomography velocity inversion method; Preferably, the multi-information constrained grid tomography velocity inversion method obtains a velocity model characterizing fine structures by iterating the depth domain velocity model, and specifically comprises the following steps: 11) Use pre-stack time migration velocity to perform time-to-depth conversion, and use VSP velocity or logging velocity and interpreted layer information to constrain the construction of the initial depth domain velocity model, complete the first round of pre-stack depth migration, output the first round of depth domain imaging volume and depth domain imaging gathers, and the migration grid of the main survey line and the contact survey line is less than or equal to 10 bins; 12) Complete the residual delay spectrum calculation and seed point picking of the first round of depth domain imaging gathers, and the picking density in the depth direction is less than or equal to 40 times the depth step; and complete the three-dimensional grid tomography velocity inversion, the spatial direction inversion grid is less than or equal to 20 facets, and the depth direction inversion grid is less than or equal to 40 times the depth step. The three-dimensional grid tomography inversion velocity model in this step is added to the initial depth domain model in step 11) to obtain a velocity model that characterizes large-scale structural features; 13) Complete the second round of pre-stack depth migration, output the second round of depth domain imaging volume and depth domain imaging gathers, and the migration grid of the main survey line and the contact survey line is less than or equal to 5 bins; 14) Complete the residual delay spectrum calculation and seed point picking of the second round of depth domain imaging gathers, and the picking density in the depth direction is less than or equal to 20 times the depth step; and complete the three-dimensional grid tomography velocity inversion, the spatial direction inversion grid is less than or equal to 10 facets, and the depth direction inversion grid is less than or equal to 20 times the depth step. The three-dimensional grid tomography inversion velocity model in this step is added to the velocity model characterizing the large-scale structural features in step 12), thereby obtaining a complex structural velocity model characterizing the medium and large scales; 15) Complete the third round of pre-stack depth migration, output the third round of depth domain imaging volume and depth domain imaging gathers, and the migration grid of the main survey line and the contact survey line is less than or equal to 2 bins; 16) Complete the calculation of the residual delay spectrum and seed point picking of 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 complete the three-dimensional grid tomographic velocity inversion, with the inversion grid in the spatial direction being less than or equal to 10 facets, and the inversion grid in the depth direction 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 structure in step 14) to obtain a velocity model representing the fine structure; 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; the iteration of the depth domain velocity model ends; 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 characterization; 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 matching the velocity model with 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.

4. 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; The type of seismic wave equation adopted in step S2, and the method for adjusting the phase of the input seismic wavelet for generating the source wave field, comprises the following steps: 21) If the wave group characteristics of the reverse time imaging result are to be consistent with those of the Kirchhoff depth imaging result, it is necessary to perform a 90-degree phase rotation on the input seismic wavelet that generates the source wave field; 22) If the reverse time imaging result is to have a wave group characteristic of approximately zero phase, the input seismic wavelet that generates the source wave field is not phase rotated; 23) 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; Preferably, the types of seismic wave equations used include first-order seismic wave equations and second-order seismic wave equations; Taking the three-dimensional acoustic wave equation as an example, the first-order seismic wave equation is: 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; The second-order seismic wave field equation is:

5. A phase-controllable and efficient reverse time migration imaging method according to claim 1 or 4, characterized in that: The specific process of completing the forward extension of the source wave field of the optimized time for any shot seismic data in the shot gather data of step S3 is as follows: In step S2, the seismic wave equation is added with the seismic wavelet generating the source wave field at x0, y0, z0. 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, and the source wave field at time 0-T is first compressed and then stored on the local disk.

6. A phase-controllable and efficient reverse time migration imaging method according to claim 1 or 4, characterized in that: 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: In this step, the actual detection point position x, y, z1(x, y) of the seismic wave equation is added to the shot data of step S1. Within the imaging space of the shot data, the calculation formula of the seismic wave equation in this step is used to perform reverse time extension of the detection point wave field of the total trace length of the shot data along the direction from large to small of the recording time of the shot data of step S1, that is, starting from the maximum trace length moment, reverse time extension has been performed to 0 moment; ① For the first-order seismic wave field equation: Where: R(x,y,z) is the wave field of the detection point; x,y,z are spatial coordinates; t represents time, and its maximum time 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 detector at the position x,y,z1(x,y), where z1(x,y) represents the depth of the detector that changes with the spatial (x,y) coordinates, including the changes in the detection point position 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; ② For the second-order seismic wave field equation: 。 7. The phase-controllable and efficient reverse time migration imaging method according to claim 5, characterized in that: The optimized source wave field propagation time length T 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; 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 the depth range.

8. 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, 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, and the specific process of completing the reverse time migration imaging processing of the shot data in the shot gather data is as follows: 31) The wave field of the detection point is reversely extended from the time corresponding to the maximum trace length of the shot data recording time in step S1. When the time Tt of the reverse time extension of the 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 the 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 the time T; 32) Then, when the time of the wave field extension of the detection point reaches the time Tt, the source wave field at the 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 the time Tt, thereby completing the reverse time imaging results of all the 0-T times; and at the same time step, the absolute value of the source wave field value is calculated to obtain the amplitude energy result at the 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 data; accumulate the amplitude energy results of all 0-T moments in step 32) to obtain the amplitude energy data body of the shot data.

9. 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, and the optimization method includes: 1) Perform amplitude compensation processing; 2) Perform low wave number noise attenuation processing; 3) According to the application requirements of 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 polarity of the output seismic wavelet after imaging due to the change of the phase of the input source wavelet in step S4, and obtain the reverse time imaging data volume processed by the final amplitude compensation; Preferably, the amplitude compensation processing method is obtained by dividing the reverse time migration data volume of the entire working 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 counted 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; Preferably, the method of setting a threshold value to limit the numerical range of the amplitude energy data volume includes: A value between [80%, 98%] of the maximum absolute value of the amplitude energy data body is used as the maximum value selection value of the modified amplitude energy data body, that is, a value in the amplitude energy data body that is greater than the maximum value selection value is set as the maximum value; A value between [1%, 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; The above processing can eliminate the abnormal processing result caused by the denominator being zero or the value being too large or too small when dividing the whole-area reverse time migration data body and the whole-area amplitude energy data body.

10. A phase-controllable and efficient reverse time migration imaging device, characterized in that: The device realizes phase-controllable and efficient reverse time migration imaging by utilizing the method described in any one of claims 1 to 9, including: The first acquisition unit is used to obtain the shot gather data after fidelity processing and the high-precision depth domain velocity model of the seismic work area; An adjustment unit, for adjusting the phase of an input seismic wavelet for generating a source wave field 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 detection point wave field with long total trace time for any shot seismic data in the shot gather data; An accumulation unit is used for applying imaging conditions at any imaging point position in the imaging space of the shot data at any time step of the optimization time in step S3, and accumulating the seismic wave field after the imaging conditions are applied within the 0-T time to complete the reverse time migration imaging processing of the shot data in the shot gather data; A second acquisition unit, based on the reverse time migration imaging processing of any single shot data in the completed shot gather data, obtains the single shot reverse time imaging data volume and the amplitude energy data volume; The first processing unit performs reverse time migration imaging processing on any single shot data in the shot gather data to obtain a single shot reverse time imaging data volume and an amplitude energy data volume; 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 whole 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