A broadband seismic wave prestack reverse time migration imaging method and device

Abstract

The present disclosure relates to a wideband seismic wave prestack reverse time migration imaging method and apparatus, including: acquiring original prestack shot gather data, depth-domain velocity model, and source data of a study area; processing the original prestack shot gather data to obtain low-frequency-preserved and fidelity-preserved wideband prestack shot gather data; using the low-frequency-preserved and fidelity-preserved wideband prestack shot gather data, depth-domain velocity model, and source data as input data, and performing reverse time migration by using an optimized reverse time migration algorithm to obtain a reverse time migration data volume; extracting an imaging gather from the reverse time migration data volume; processing the imaging gather to obtain a reverse time migration imaging result. This is to solve the problems that in previous reverse time migration imaging, all wavefield frequency information of the shot gather data was not utilized, there was no optimized consideration of the frequency bandwidth in the migration parameters during migration, the migration frequency was limited by a strict numerical dispersion relationship, the resolution of the imaging result was not high, and the ability to finely depict seismic complex geological targets was limited.

Description

Technical Field

[0001] The present disclosure relates to the technical field of seismic exploration data processing, and particularly to a broadband seismic wave prestack reverse time migration imaging method and apparatus. Background Art

[0002] With the research and development and large-scale application of seismic wave reverse time imaging technology, significant geological effects have been achieved in accurate imaging of complex wave fields and complex structures, and it has received extensive attention and hot research in the international geophysical community. The reverse time imaging technology is a technology that uses the two-way seismic wave wave equation and combines high-order approximation numerical discretization algorithms, which can effectively solve key technical problems such as multi-path propagation, multiple wave migration, imaging of large and steep dips, imaging of turning waves, imaging of prism waves, and amplitude-preserved seismic imaging. It is suitable for high-precision imaging of complex wave fields with relatively drastic velocity and structural changes, and effectively solves the problem that Kirchhoff integral method depth migration and one-way wave equation depth migration cannot achieve accurate migration of complex wave fields. Therefore, it has obvious advantages in seismic imaging technology. With the rapid development of high-performance computing technologies such as CPU / GPU and the large-scale application of supercomputing platforms and mass parallel storage technologies in recent years, the reverse time imaging technology has broken through a series of technical bottlenecks such as mass storage and huge computational volume, realized large-scale production applications, and effectively supported the exploration and development of complex geological targets.

[0003] Existing problems: The industry generally believes that reverse time migration can only solve the problem of seismic imaging of complex structures. Strictly according to the numerical dispersion relationship, under the existing acquisition bin and seismic data conditions, the maximum migration frequency of reverse time migration is usually selected as 20 - 30 Hz, while the actual seismic data frequency band is 5 - 100 Hz. Therefore, all wave field frequency information of prestack seismic data cannot be fully utilized. Using a higher migration frequency will introduce strong numerical dispersion interference, thus greatly reducing the signal-to-noise ratio of the imaging result. To improve the frequency band and resolution of the reverse time migration imaging result, it is necessary to pay the price of a more huge computational volume and cost. Therefore, how to efficiently and low-costly implement broadband reverse time migration processing and improve the ability to depict details of complex geological targets is the main difficult problem faced by current high-precision reverse time migration. Summary of the Invention

[0004] The present disclosure provides a broadband seismic wave prestack reverse time migration imaging method and apparatus to solve the problems that when performing reverse time migration imaging in the past, all wave field frequency information of shot gather data was not utilized, there was no optimized consideration of the frequency bandwidth in the migration parameters during migration, the migration frequency was limited by the strict numerical dispersion relationship, the resolution of the imaging result was not high, and the ability to finely depict the seismic conditions of complex geological targets was limited.

[0005] According to one aspect of the present disclosure, there is provided a broadband seismic wave prestack reverse time migration imaging method, including:

[0006] Obtain the original prestack shot gather data, depth-domain velocity model, and source data in the study area;

[0007] Perform denoising, spherical spreading compensation, surface-consistent energy compensation, and deconvolution on the original prestack shot gather data to obtain low-frequency-preserved, fidelity-preserved, wide-band prestack shot gather data;

[0008] Utilize the depth-domain velocity model to optimize the parameters of the reverse time migration algorithm, obtaining the reverse time migration algorithm with optimized parameters;

[0009] Use the low-frequency-preserved, fidelity-preserved, wide-band prestack shot gather data, depth-domain velocity model, and source data as input data, and perform reverse time migration using the reverse time migration algorithm with optimized parameters to obtain a reverse time migration data volume;

[0010] Extract the imaging gather from the reverse time migration data volume;

[0011] Perform truncation stacking and low wavenumber noise removal processing on the imaging gather to obtain the reverse time migration imaging result.

[0012] Preferably, before utilizing the depth-domain velocity model to optimize the parameters of the reverse time migration algorithm, it further includes:

[0013] Perform prestack interpolation processing on the low-frequency-preserved, fidelity-preserved, wide-band prestack shot gather data in the common shot point domain, and perform anti-spherical spreading compensation processing on the prestack shot gather data after the interpolation processing.

[0014] Preferably, before obtaining the depth-domain velocity model, it further includes: establishing a depth-domain velocity model;

[0015] The method for establishing the depth-domain velocity model includes:

[0016] Obtain an isotropic velocity model vel or an anisotropic parameter model;

[0017] Perform grid tomography global inversion on the isotropic velocity model vel or the anisotropic parameter model to obtain a depth-domain velocity model;

[0018] Wherein, the depth step of the grid tomography global inversion is within a first threshold, and the inversion grid step in the X and Y directions is less than a first predetermined value and greater than 2 bin sizes.

[0019] Preferably, the method for utilizing the depth-domain velocity model to optimize the parameters of the reverse time migration algorithm to obtain the reverse time migration algorithm with optimized parameters includes:

[0020] Set the time difference accuracy of the reverse time migration numerical discretization algorithm to be greater than or equal to a first predetermined value, and the space difference accuracy to be greater than or equal to a second predetermined value;

[0021] Under the differential accuracy condition, using the depth-domain velocity model, perform inverse time migration parameter optimization to obtain an inverse time migration algorithm with optimized parameters.

[0022] Preferably, the method for performing inverse time migration parameter optimization using the depth-domain velocity model under the differential accuracy condition includes:

[0023] Obtain the minimum velocity in the depth-domain velocity model;

[0024] Determine the number of grid nodes occupied by each seismic wave wavelength under the differential accuracy;

[0025] Set the maximum bin and minimum depth step of the migration;

[0026] According to the number of grid nodes, the minimum velocity, the maximum bin, and the minimum depth step, determine the maximum migration frequency allowed for reverse time migration in the spatial direction and the maximum migration frequency allowed for reverse time migration in the depth direction;

[0027] According to the maximum migration frequency allowed for reverse time migration in the spatial direction and the maximum migration frequency allowed for reverse time migration in the depth direction, determine the maximum frequency parameter used for wide-band reverse time migration.

[0028] Preferably, the method for using the optimized inverse time migration algorithm to perform inverse time migration on the low-frequency-preserved, fidelity-preserved, wide-band pre-stack shot gather data, depth-domain velocity model, and source data as input data to obtain an inverse time migration data volume includes:

[0029] Using the optimized inverse time migration algorithm, in the depth-domain velocity model, perform wavefield continuation on the low-frequency-preserved, fidelity-preserved, wide-band pre-stack shot gather data and source data to generate the wavefield of the receiving points of the pre-stack shot gather data and the source wavefield, and perform up-going wave and down-going wavefield separation processing on the wavefield of the receiving points of the pre-stack shot gather data and the source wavefield at each moment;

[0030] According to the separated wavefields, determine the inverse time imaging condition at each time step for each imaging position;

[0031] Accumulate the inverse time imaging conditions at all time steps corresponding to each imaging position to obtain the imaging point value at this imaging position, and all the imaging point values constitute the inverse time migration data volume.

[0032] Preferably, the method for determining the inverse time imaging condition at each time step for each imaging position according to the separated wavefields includes:

[0033] According to the separated wavefields, use the inverse time imaging condition determination formula to determine the inverse time imaging condition at each time step for each position in the imaging space;

[0034] Among them, the formula for determining the inverse-time imaging condition includes:

[0035] I = R 上 *S 上 +R 下 *S 下 ;

[0036] In the formula: R 上 is the up-going wave of the wave field at the detection point, S 上 is the up-going wave of the wave field at the source, R 下 is the down-going wave of the wave field at the detection point, S 下 is the down-going wave of the wave field at the source.

[0037] Preferably, the method for extracting the imaging gather in the inverse-time migration data volume includes:

[0038] In the inverse-time migration data volume, extract the data information of the trace headers of each shot to form a common image point super gather;

[0039] According to the set number of azimuths, starting offset, and offset interval, reconstruct the data traces of the imaging point super gather to obtain the imaging gather.

[0040] Preferably, the method for muting stacking and removing low wavenumber noise processing includes:

[0041] Pick up the muting function on the imaging gather at any position in the fully covered area of the study area, and apply the muting function to the imaging gather to obtain the imaging gather after muting the stretched distorted wave field;

[0042] Accumulate the values of the imaging gather after muting the stretched distorted wave field along the same depth direction to obtain the inverse-time migration stacked data volume;

[0043] Use the diffusion filtering method to perform relative amplitude-preserving suppression processing on the low wavenumber noise of the inverse-time migration stacked data volume to obtain the inverse-time migration imaging result.

[0044] According to one aspect of the present disclosure, there is provided a wideband seismic wave pre-stack inverse-time migration imaging device, including:

[0045] An acquisition unit for acquiring the original pre-stack shot gather data, velocity model in the depth domain, and source data of the study area;

[0046] A pre-stack shot gather data processing unit for denoising, spherical spreading compensation, surface consistent energy compensation, and deconvolution processing on the original pre-stack shot gather data to obtain low-frequency-preserved and fidelity wideband pre-stack shot gather data;

[0047] A parameter optimization unit, configured to optimize the parameters of the reverse time migration algorithm by using the depth-domain velocity model, so as to obtain a reverse time migration algorithm with optimized parameters;

[0048] A reverse time migration unit, configured to use the low-frequency-preserved and fidelity-wideband pre-stack gather data, the depth-domain velocity model, and the source data as input data, and perform reverse time migration by using the reverse time migration algorithm with optimized parameters to obtain a reverse time migration data volume;

[0049] An imaging gather extraction unit, configured to extract an imaging gather from the reverse time migration data volume;

[0050] An imaging result generation unit, configured to perform excision stacking and low wavenumber noise removal processing on the imaging gather to obtain a reverse time migration imaging result.

[0051] The present invention has at least the following beneficial effects:

[0052] The present disclosure provides a wideband seismic wave pre-stack reverse time migration imaging method and apparatus. By using the obtained wideband low-frequency-preserved and fidelity-wideband pre-stack gather data of the study area and the refined velocity model as the input data for reverse time migration, all frequency band information of the pre-stack data is fully utilized. Then, combined with the reverse time migration parameter optimization method based on the numerical dispersion relationship, wideband reverse time migration processing is realized, and excision and optimized stacking processing are performed on the migrated gathers. After comprehensive application, the frequency band is increased by more than 1 time compared with the conventional method. Without increasing the computational complexity and cost, the fine characterization ability of complex geological targets is finally improved efficiently and with high precision. Description of the Drawings

[0053] The drawings herein are incorporated into the specification and form a part of this specification. These drawings show embodiments consistent with the present disclosure and, together with the specification, are used to explain the technical solutions of the present disclosure.

[0054] Figure 1 A flowchart showing a wideband seismic wave pre-stack reverse time migration imaging method according to an embodiment of the present disclosure.

[0055] Figure 2 Showing a conventional reverse time migration imaging result of the international standard BP theoretical model strictly in accordance with the numerical dispersion relationship according to an embodiment of the present disclosure.

[0056] Figure 3 Showing a wideband reverse time migration imaging result of the international standard two-dimensional BP theoretical model obtained by using the wideband seismic wave pre-stack reverse time migration imaging method according to an embodiment of the present disclosure.

[0057] Figure 4 Showing a reverse time migration imaging result obtained strictly in accordance with the numerical dispersion relationship for seismic data in a certain study area in the middle and shallow layers of the Songliao Basin according to an embodiment of the present disclosure.

[0058] Figure 5 Shows the broadband reverse time migration imaging result obtained by using the broadband seismic wave pre-stack reverse time migration imaging method for seismic data in a certain study area in the middle and shallow layers of the Songliao Basin according to an embodiment of the present disclosure. Detailed implementation manners

[0059] Various exemplary embodiments, features and aspects of the present disclosure will be described in detail below with reference to the accompanying drawings. The same reference numerals in the drawings denote elements having the same or similar functions. Although various aspects of the embodiments are shown in the drawings, the drawings are not necessarily drawn to scale unless otherwise specified.

[0060] The term "exemplary" used herein means "serving as an example, embodiment or illustration". Any embodiment described as "exemplary" herein is not necessarily to be construed as superior or better than other embodiments.

[0061] The term "and / or" in this article is merely a description of the associated relationship of the associated objects, indicating that there can be three relationships. For example, A and / or B can represent: A exists alone, A and B exist simultaneously, and B exists alone. In addition, the term "at least one" in this article means any one of a plurality or any combination of at least two of a plurality. For example, including at least one of A, B, and C can represent including any one or more elements selected from the set composed of A, B, and C.

[0062] In addition, for better illustration of the present disclosure, numerous specific details are given in the following detailed implementation manners. Those skilled in the art should understand that the present disclosure can also be implemented without some specific details. In some instances, methods, means, elements and circuits well known to those skilled in the art are not described in detail so as to highlight the gist of the present disclosure.

[0063] Figure 1 Shows the flowchart of the broadband seismic wave pre-stack reverse time migration imaging method according to an embodiment of the present disclosure; Figure 2 Shows the conventional reverse time migration imaging result of the international standard BP theoretical model in strict accordance with the numerical dispersion relationship according to an embodiment of the present disclosure; Figure 3 Shows the broadband reverse time migration imaging result of the international standard two-dimensional BP theoretical model obtained by using the broadband seismic wave pre-stack reverse time migration imaging method according to an embodiment of the present disclosure; Figure 4 Shows the reverse time migration imaging result obtained in strict accordance with the numerical dispersion relationship for the seismic data in a certain study area in the middle and shallow layers of the Songliao Basin according to an embodiment of the present disclosure; Figure 5 Shows the broadband reverse time migration imaging result obtained by using the broadband seismic wave pre-stack reverse time migration imaging method for the seismic data in a certain study area in the middle and shallow layers of the Songliao Basin according to an embodiment of the present disclosure. As Figures 1-5As shown in the figure, a wideband seismic wave pre-stack reverse time migration imaging method includes the following steps: Step S01: Obtain the original pre-stack shot gather data, depth-domain velocity model, and source data of the study area; Step S02: Denoise, spherical spreading compensation, surface consistent energy compensation, and deconvolution processing are performed on the original pre-stack shot gather data to obtain a low-frequency and high-fidelity wideband pre-stack shot gather data; Step S03: Using the depth-domain velocity model, optimize the parameters of the reverse time migration algorithm to obtain the reverse time migration algorithm with optimized parameters; Step S04: Use the low-frequency and high-fidelity wideband pre-stack shot gather data, depth-domain velocity model, and source data as input data, and perform reverse time migration using the reverse time migration algorithm with optimized parameters to obtain a reverse time migration data volume; Step S05: Extract the imaging gather from the reverse time migration data volume; Step S06: Perform excision stacking and low wavenumber noise removal processing on the imaging gather to obtain the reverse time migration imaging result.

[0064] The wideband seismic wave pre-stack reverse time migration imaging method provided by the embodiment of the present invention specifically includes the following steps:

[0065] Step S01: Obtain the original pre-stack shot gather data, depth-domain velocity model, and source data of the study area.

[0066] Step S02: Denoise, spherical spreading compensation, surface consistent energy compensation, and deconvolution processing are performed on the original pre-stack shot gather data to obtain a low-frequency and high-fidelity wideband pre-stack shot gather data.

[0067] In the embodiment of the present disclosure, the denoising process includes pre-stack fidelity denoising such as surface wave noise removal, and a suitable denoising method is selected according to the specific situation of the geological data in the study area; the spherical spreading compensation process and the surface consistent energy compensation process are used to compensate for the energy absorbed and attenuated by the formation; the deconvolution process includes a combination of surface consistent deconvolution and predictive deconvolution, which is used to improve the resolution of the pre-stack shot gather data and broaden the frequency band of the pre-stack shot gather data.

[0068] Step S03: Using the depth-domain velocity model, optimize the parameters of the reverse time migration algorithm to obtain the reverse time migration algorithm with optimized parameters.

[0069] In the present disclosure, before using the depth-domain velocity model to optimize the parameters of the reverse time migration algorithm, it further includes: performing pre-stack interpolation processing on the low-frequency and high-fidelity wideband pre-stack shot gather data in the common shot point domain, and performing anti-spherical spreading compensation processing on the pre-stack shot gather data after the interpolation processing.

[0070] In the embodiments of the present disclosure, the interpolation process includes NMO processing, pre-stack interpolation processing, and de-NMO processing; the interpolation process is used to reduce the geophone point spacing and geophone line spacing within any shot data, thereby increasing the fold number of the gather data, increasing the trace density and data volume of each shot, and meeting the requirements of high-precision reverse-time imaging for small bins, without the problem of migration arc interference introduced when imaging a smaller bin with the originally acquired bin shot gather data.

[0071] Perform spherical divergence compensation processing on the pre-stack shot gather data after interpolation processing, and finally obtain the pre-stack shot gather data Shotgather for input data of wide-band reverse-time migration processing.

[0072] In the present disclosure, before obtaining the depth-domain velocity model, it further includes: establishing a depth-domain velocity model; the method for establishing the depth-domain velocity model includes: obtaining an isotropic velocity model vel or an anisotropic parameter model; performing global grid tomography inversion on the isotropic velocity model vel or the anisotropic parameter model to obtain a depth-domain velocity model; wherein, the depth step of the global grid tomography inversion is within a first threshold, and the inversion grid step sizes in the X and Y directions are less than a first predetermined value and greater than twice the bin size.

[0073] In the embodiments of the present disclosure, the velocity model can be an isotropic velocity model vel or an anisotropic parameter model; the velocity model includes one or several of anisotropic velocity veldata, anisotropic parameter field deltadata, anisotropic parameter field epsdata, dip angle field dipxdata in the x direction, and dip angle field dipydata in the y direction.

[0074] Reduce the inversion grid size of the isotropic velocity model vel or the anisotropic parameter model through the global grid tomography inversion method. Among them, the first threshold is 200m - 5m, preferably 50m - 5m; the value range of the first predetermined value is 400m - 100m, preferably 200m; the depth step of the global grid tomography inversion is less than 50m and greater than 5m, and the inversion grid step sizes in the X and Y directions are less than 400m and greater than twice the bin size. Finally, obtain a depth-domain velocity model for parameter optimization of wide-band reverse-time migration processing and as its input data.

[0075] Step S03: Use the depth-domain velocity model to optimize the parameters of the reverse-time migration algorithm, and obtain the reverse-time migration algorithm with optimized parameters.

[0076] In the present disclosure, the method for optimizing the parameters of the reverse time migration algorithm by using a depth-domain velocity model to obtain an optimized reverse time migration algorithm includes: setting the time difference accuracy of the numerical discretization algorithm of reverse time migration to be greater than or equal to a first predetermined value, and the spatial difference accuracy to be greater than or equal to a second predetermined value; under the condition of the difference accuracy, using the depth-domain velocity model to perform reverse time migration parameter optimization to obtain an optimized reverse time migration algorithm.

[0077] In an embodiment of the present disclosure, the first predetermined value is 2nd order, and the second predetermined value is 8th order. The time difference accuracy of the numerical discretization algorithm of reverse time migration is higher than 2nd order, and the spatial difference accuracy is higher than 8th order.

[0078] In the present disclosure, the method for performing reverse time migration parameter optimization by using a depth-domain velocity model under the condition of the difference accuracy includes: obtaining the minimum velocity velmin in the depth-domain velocity model; determining the number of grid nodes occupied by each seismic wave wavelength under the difference accuracy; setting the maximum bin binmax and the minimum depth step depthmin of the migration; according to the number of grid nodes, the minimum velocity velmin, the maximum bin binmax, and the minimum depth step depthmin, determining the maximum migration frequency allowed for reverse time migration in the spatial direction and the maximum migration frequency allowed for reverse time migration in the depth direction; according to the maximum migration frequency allowed for reverse time migration in the spatial direction and the maximum migration frequency allowed for reverse time migration in the depth direction, determining the maximum frequency parameter fmax adopted for wideband reverse time migration.

[0079] In an embodiment of the present disclosure, the number of grid nodes occupied by the wavelength of each seismic wave is a quantity closely related to the difference accuracy. The higher the difference accuracy, the smaller the number of nodes required for each wavelength can be, but it must be greater than 2, that is, to satisfy the sampling theorem, at least 2 points are required for each wavelength to reconstruct or recover a complete seismic waveform; the more points occupied by each wavelength, such as 10 or higher, the more accurate the reconstructed or recovered seismic waveform.

[0080] In this embodiment, the selected time difference accuracy is 2nd order, and the spatial difference accuracy is 8th order. Under this condition, it is at least necessary to ensure that the number of grid nodes occupied by each wavelength is greater than 2.5.

[0081] If the number of grid nodes is selected as 2.5, according to the minimum velocity velmin of the depth-domain velocity model, the set maximum bin binmax and the minimum depth step depthmin of the migration, according to formula (1) and formula (2), determine the maximum migration frequency fmax1 allowed for reverse time migration in the spatial direction and the maximum migration frequency fmax2 allowed for reverse time migration in the depth direction.

[0082] fmax1 = velmin / (2.5 * binmax); (1)

[0083] fmax2 = velmin / (2.5 * depthmin); (2)

[0084] Under normal circumstances, fmax2 > fmax1 and binmax > depthmin. Then the maximum frequency parameter in the broadband reverse time migration parameter is: fmax = 2 * fmax1. When the calculated result fmax >= fmax2, then fmax = fmax2.

[0085] Step S04: Use the low-frequency-preserving and fidelity-preserving broadband pre-stack gather data, the depth-domain velocity model, and the source data as input data, and perform reverse time migration using the optimized reverse time migration algorithm to obtain a reverse time migration data volume.

[0086] In the present disclosure, the method of using the low-frequency-preserving and fidelity-preserving broadband pre-stack gather data, the depth-domain velocity model, and the source data as input data, and performing reverse time migration using the optimized reverse time migration algorithm to obtain a reverse time migration data volume includes: using the optimized reverse time migration algorithm, in the depth-domain velocity model, performing wavefield continuation on the low-frequency-preserving and fidelity-preserving broadband pre-stack gather data and the source data to generate the wavefield of the receiving points of the pre-stack gather data and the source wavefield, and performing up-going wave and down-going wave field separation processing on the wavefield of the receiving points of the pre-stack gather data and the source wavefield at each moment; determining the reverse time imaging condition at each time step for each imaging position according to the separated wavefield; accumulating the reverse time imaging conditions at all time steps corresponding to each imaging position to obtain the imaging point value at this imaging position, and all the imaging point values constitute the reverse time migration data volume.

[0087] In the embodiment of the present disclosure, the method of performing wavefield continuation on the low-frequency-preserving and fidelity-preserving broadband pre-stack gather data to generate the wavefield of the receiving points includes: using the maximum frequency parameter fmax in the optimized reverse time migration algorithm obtained in step S03 to perform low-pass filtering on the low-frequency-preserving and fidelity-preserving broadband pre-stack gather data to obtain the pre-stack gather data after low-pass filtering; using the discretized seismic wave equation in the reverse time migration algorithm to perform reverse time wavefield continuation (propagating downward) on the pre-stack gather data after low-pass filtering in the depth-domain velocity model to generate the wavefield of the receiving points.

[0088] If the adopted frequency is greater than fmax, then the frequency information greater than fmax in the gather data will perform wavefield continuation in the form of interference, resulting in numerical dispersion artifacts, thereby interfering with the accuracy of the final wavefield of the receiving points. Therefore, before performing wavefield continuation, it is necessary to perform low-pass filtering using fmax first.

[0089] In an embodiment of the present disclosure, a method for wavefield continuation of source wave data to generate a source wavefield includes: obtaining a band-limited wavelet as a source according to the maximum frequency parameter fmax and the band-limited wavelet function formula; and using the discretized seismic wave equation in the reverse time migration algorithm to generate a source wavefield through forward wavefield continuation (propagating downward) in the depth-domain velocity model.

[0090] Among them, the band-limited wavelet function formula includes:

[0091]

[0092] In the formula, f H is the upper limit frequency of the band-limited wavelet, that is, fmax, and f L is the lower limit frequency of the band-limited wavelet, and f L is not equal to 0, and t is time.

[0093] If the adopted frequency is greater than fmax, then the frequency information greater than fmax in the source wavefield will be wavefield continued in the form of interference, resulting in a numerical dispersion artifact, thus interfering with the accuracy of the final source wavefield.

[0094] In an embodiment of the present disclosure, during the imaging process of each shot, the forward-propagating wavefield of the source is separated into an up-going wave S 上 and a down-going wave S 下 , and the reverse-time continued wavefield of the geophone is separated into an up-going wave R 上 and a down-going wave R down.

[0095] In the present disclosure, the method for determining the reverse-time imaging condition at each time step for each imaging position according to the separated wavefields includes: determining the reverse-time imaging condition at each time step for each position in the imaging space according to the separated wavefields by using the reverse-time imaging condition determination formula; among them, the reverse-time imaging condition determination formula includes:

[0096] I = R 上 * S 上 + R 下 * S 下 ;

[0097] In the formula: R 上 is the up-going wave of the geophone wavefield, S 上 is the up-going wave of the source wavefield, R 下 is the down-going wave of the geophone wavefield, and S 下 is the down-going wave of the source wavefield.

[0098] In an embodiment of the present disclosure, the reverse-time imaging condition determination formula is applied to process at each time step for each imaging position in the imaging space to obtain the reverse-time imaging condition I corresponding to each time step.

[0099] In the embodiments of the present disclosure, after obtaining the inverse-time imaging condition I for each time step corresponding to each imaging position in the imaging space, an accumulation process is performed on the inverse-time imaging conditions corresponding to all time steps at each imaging position in the imaging space to obtain the imaging point value corresponding to this imaging position. All the imaging point values constitute the inverse-time migration data volume, and thus the inverse-time migration of all low-frequency-preserved, fidelity-preserved, wide-band, prestack shot gather data in the study area is completed.

[0100] Step S05: Extract an imaging gather from the inverse-time migration data volume.

[0101] In the present disclosure, the method for extracting an imaging gather from the inverse-time migration data volume includes:

[0102] In the inverse-time migration data volume, extract the data information of the trace headers for each shot to form a common image point super gather; according to the set number of azimuths, starting offset, and offset interval, perform data trace reconstruction on the common image point super gather to obtain the imaging gather.

[0103] In the embodiments of the present disclosure, quickly read all the data belonging to this imaging position from the inverse-time migration data volume of all the shot gathers in the study area obtained in step S04 to form a common image point super gather at this spatial position. Then, perform data trace reconstruction according to the preset number of azimuths, starting offset, and offset interval.

[0104] Among them, data trace reconstruction includes: accumulating the data traces belonging to a certain azimuth angle and within a certain offset interval in this azimuth angle. If there are only 0 data traces within this offset interval range, zero-fill the data at this position. If there are more than 1 data traces within this offset interval range, accumulate the sampled point values along the same depth direction, thereby obtaining the offset imaging gathers at all imaging positions in the study area.

[0105] Among them, the extracted imaging gather can be a single-azimuth offset imaging gather or an omnidirectional offset imaging gather. The number of azimuth angles of the imaging gather can be 6 azimuths, and the offset interval can be 50 - 100 m, which can be specifically determined according to the geological requirements of the actual study area.

[0106] Step S06: Perform truncation stacking and low wavenumber noise removal processing on the imaging gather to obtain the inverse-time migration imaging result.

[0107] In the present disclosure, the method for resection superposition and low wavenumber noise removal processing includes: picking up a resection function on the imaging gather at any position in the full-coverage position of the study area, applying the resection function to the imaging gather to obtain an imaging gather after resection of the stretched distorted wavefield; accumulating the values of the imaging gather after resection of the stretched distorted wavefield along the same depth direction to obtain an inverse time migration superposition data volume; performing inverse time migration low wavenumber noise relative amplitude-preserving suppression processing on the inverse time migration superposition data volume by using a diffusion filtering method to obtain an inverse time migration imaging result.

[0108] In an embodiment of the present disclosure, the specific formula of the diffusion filtering method includes:

[0109] U 0 = StackRegu, t = 0, 1, 2,..., N;

[0110] In the formula: α is the iterative filtering coefficient, N is the total number of iterative filtering times, t represents the current iterative time, U 0 is the initial migration superposition data volume, U t+1 is the migration superposition data volume after the (t + 1)-th iteration.

[0111] In an embodiment of the present disclosure, taking the international standard two-dimensional BP theoretical model as an example, original pre-stack shot gather data, a depth-domain velocity model, and source data are acquired. Denoising processing, spherical spreading compensation processing, surface-consistent energy compensation processing, and deconvolution processing are performed on the original pre-stack shot gather data to obtain low-frequency-preserving, fidelity-preserving, wide-frequency pre-stack shot gather data.

[0112] Using the depth-domain velocity model, the parameters of the inverse time migration algorithm are optimized. Among them, the time difference accuracy of the numerical discretization algorithm for inverse time migration is set to 2nd order, and the spatial difference accuracy is set to 16th order accuracy; the minimum velocity in the velocity model is obtained as 1875 m / s, the maximum bin for migration is selected as 25 m, the minimum depth step is 12.5 m, and the number of grid nodes occupied by each wavelength is ensured to be 2.5 points or more under the set difference accuracy conditions. According to formulas (1) and (2), the maximum migration frequency allowed for inverse time migration in the spatial direction is determined to be 30 Hz, and the maximum migration frequency allowed for inverse time migration in the depth direction is determined to be 60 Hz. Thus, the maximum frequency parameter fmax used for wide-frequency inverse time migration is obtained as 60 Hz.

[0113] Using the reverse time migration algorithm with optimized parameters, in the depth domain velocity model, wavefield continuation is performed on the low-frequency-preserved, fidelity-preserved, wide-band pre-stack shot gather data and source data to generate the wavefield of the receiver points and the source wavefield of the pre-stack shot gather data. The wavefield at each moment is processed for separating the up-going wave and down-going wave fields. According to the separated wave fields, the reverse time imaging condition using selective wavefield components for correlation imaging is adopted to eliminate the low wavenumber background noise in the imaging process, and a reverse time migration data volume is obtained.

[0114] An imaging gather is extracted from the reverse time migration data volume. This imaging gather can be a single azimuth offset imaging gather with an offset interval of 50 m. The excision function is picked up and applied to the offset imaging gathers at all positions in the study area to obtain the imaging gathers after excision of the stretched distorted wavefield. The values of the excised imaging gathers along the same depth direction are accumulated to obtain the reverse time migration stacked data volume. The diffusion filtering method is used for relatively amplitude-preserving suppression of the low wavenumber noise in reverse time migration to obtain the final wide-band reverse time migration imaging result.

[0115] Figure 2 For the conventional reverse time migration imaging result of the international standard BP theoretical model strictly in accordance with the numerical dispersion relationship, where the maximum bin of the migration grid is 25 m, the depth step is 12.5 m, and the maximum migration frequency is 24 Hz; Figure 3 For the wide-band reverse time migration imaging result of the international standard two-dimensional BP theoretical model according to the wide-band seismic wave pre-stack reverse time migration imaging method of the present disclosure, where the maximum bin of the migration grid is 25 m, the depth step is 12.5 m, and the maximum migration frequency is 60 Hz. From Figure 2 、 Figure 3 By comparison, it can be seen that both sets of migration parameters can achieve accurate wavefield imaging of complex geological targets. Among them, Figure 3 The wide-band reverse time migration imaging result makes full use of the wavefield information and realizes the fine characterization of the geological details of the complex wavefield.

[0116] In the embodiment of the present disclosure, taking the seismic data of a certain study area in the middle and shallow layers of the Songliao Basin as an example, the original pre-stack shot gather data of the study area is obtained. Through pre-stack fidelity denoising processing such as surface wave noise removal, spherical divergence compensation and surface consistent energy compensation processing, and the combined processing of surface consistent deconvolution and predictive deconvolution, the pre-stack shot gather data with broadened frequency band under the condition of low-frequency preservation and fidelity preservation is obtained.

[0117] Pre-stack interpolation processing is performed on the low-frequency-preserved, fidelity-preserved, wide-band pre-stack shot gather data in the common shot point domain, specifically including moveout correction processing, interpolation processing and inverse moveout correction processing. The interval between receiver points and receiver line intervals within each shot data are reduced by half, thereby increasing the fold number and meeting the requirements of high-precision reverse time imaging for small bins; anti-spherical divergence compensation processing is performed on the pre-stack shot gather data, and finally the low-frequency-preserved, fidelity-preserved, wide-band pre-stack shot gather data Shotgather for wide-band reverse time migration processing is obtained.

[0118] Obtain a high-precision depth-domain velocity model for the study area. This model can be an isotropic velocity model vel. The model is obtained through a global tomographic inversion method that reduces the inversion grid size. The inversion depth step is 50 m, and the inversion grid sizes in the X and Y directions are 200 m.

[0119] Use the low-frequency-preserving, fidelity-preserving, wide-band pre-stack shot gather data and the high-precision depth-domain velocity model as the input data for reverse time migration, and complete the reverse time imaging process of the data using a high-precision reverse time migration algorithm with optimized migration parameters. Among them, the time difference accuracy is 2nd order and the spatial difference accuracy is 16th order when optimizing the reverse time migration algorithm parameters; the minimum velocity in the obtained depth-domain velocity model is 1500 m / s, the maximum bin size for migration is selected as 20 m, and the minimum depth step is 5 m; the number of grid nodes occupied by each wavelength is guaranteed to be 2.5 points or more under the set difference accuracy conditions. According to formulas (1) and (2), the maximum migration frequency allowed for reverse time migration in the spatial direction is 30 Hz, and the maximum migration frequency allowed for reverse time migration in the depth direction is 60 Hz. Thus, the maximum frequency parameter fmax used for wide-band reverse time migration is 60 Hz.

[0120] Using the reverse time migration algorithm with optimized parameters, in the depth-domain velocity model, perform wavefield extrapolation on the low-frequency-preserving, fidelity-preserving, wide-band pre-stack shot gather data and the source data to generate the wavefield of the geophone points in the pre-stack shot gather data and the source wavefield; during the imaging process of each shot, separate the forward-propagating source wavefield into the up-going wave S 上 and the down-going wave S 下 , separate the reverse-time extrapolated wavefield of the geophone points into the up-going wave R 上 and the down-going wave R 下 , apply the reverse time imaging condition I = R 上 *S 上 +R 下 *S 下 for processing at each imaging position in the imaging space for each time step, and perform cumulative processing on the reverse time imaging conditions corresponding to all time steps at this imaging position in the imaging space to obtain the reverse time migration data volume.

[0121] Extract the imaging gather from the reverse time migration data volume. This imaging gather is a 6-azimuth offset imaging gather with an offset interval of 50 m. Extract the offset imaging gather at any position in the fully covered area of the work area, pick up the excision function, and apply it to the 6-azimuth offset imaging gathers in the entire work area to obtain the imaging gather after removing the stretched and distorted wavefield. Accumulate the values of the imaging gather after excision along the same depth direction to obtain the reverse time migration stacked data volume. Use the diffusion filtering method to perform relative amplitude suppression processing on the low wavenumber noise in the reverse time migration to obtain the final wide-band reverse time migration imaging result.

[0122] Figure 4 It is the reverse time migration imaging result obtained when the seismic data of a certain research area in the middle and shallow layers of the Songliao Basin determines the migration parameters strictly according to the numerical dispersion relationship. Among them, the migration grid is 20m, the depth step is 5m, and the maximum migration frequency is 30Hz; Figure 5 It is the broadband reverse time migration imaging result obtained from the seismic data of a certain research area in the middle and shallow layers of the Songliao Basin according to the broadband seismic wave pre-stack reverse time migration imaging method of the present disclosure. Among them, the migration grid is 20m, the depth step is 5m, and the maximum migration frequency is 60Hz. Comparison Figure 4 、 Figure 5 It can be seen that both sets of migration parameters can achieve accurate wave field imaging of complex geological targets in the Songliao Basin. Among them, Figure 5 the broadband reverse time migration imaging result makes full use of the broadband input seismic wave field information and realizes the fine characterization of complex geological targets.

[0123] It can be understood that the above-mentioned various method embodiments mentioned in the present disclosure can be combined with each other to form a combined embodiment without violating the principle logic. Due to space limitations, the present disclosure will not elaborate further.

[0124] The execution subject of the broadband seismic wave pre-stack reverse time migration imaging method can be a broadband seismic wave pre-stack reverse time migration imaging device. For example, the broadband seismic wave pre-stack reverse time migration imaging method can be executed by a terminal device, a server, or other processing devices. Among them, the terminal device can be a user equipment (UE), a mobile device, a user terminal, a terminal, a cellular phone, a cordless phone, a personal digital assistant (PDA), a handheld device, a computing device, a vehicle-mounted device, a wearable device, etc. In some possible implementation manners, the broadband seismic wave pre-stack reverse time migration imaging method can be implemented by a processor calling computer-readable instructions stored in a memory.

[0125] Those skilled in the art can understand that in the above method of the specific implementation manner, the writing order of each step does not mean a strict execution order and does not constitute any limitation to the implementation process. The specific execution order of each step should be determined according to its function and possible internal logic.

[0126] The present disclosure provides a wideband seismic wave pre-stack reverse time migration imaging device, comprising: an acquisition unit configured to acquire original pre-stack shot gather data, a depth-domain velocity model, and source data of a study area; a pre-stack shot gather data processing unit configured to perform denoising, spherical spreading compensation, surface-consistent energy compensation, and deconvolution processing on the original pre-stack shot gather data to obtain low-frequency-preserved and fidelity-preserved wideband pre-stack shot gather data; a parameter optimization unit configured to optimize reverse time migration algorithm parameters by using the depth-domain velocity model to obtain an optimized reverse time migration algorithm; a reverse time migration unit configured to use the low-frequency-preserved and fidelity-preserved wideband pre-stack shot gather data, the depth-domain velocity model, and the source data as input data and perform reverse time migration by using the optimized reverse time migration algorithm to obtain a reverse time migration data volume; an imaging gather extraction unit configured to extract imaging gathers from the reverse time migration data volume; and an imaging result generation unit configured to perform excision stacking and low wavenumber noise removal processing on the imaging gathers to obtain a reverse time migration imaging result.

[0127] In some embodiments, the functions or modules and units included in the device provided by the embodiments of the present disclosure can be used to execute the methods described in the method embodiments above. The specific implementation can refer to the description of the method embodiments above. For the sake of brevity, it will not be repeated here.

[0128] Aiming at the problems of the existing reverse time migration technology that do not optimize the frequency bandwidth in terms of the quality of the input shot gather data and migration parameters, strictly rely on the numerical dispersion relationship to determine the migration parameters, resulting in limited wavefield information of the seismic data used in reverse time migration, reduced imaging resolution, large increase in computational amount, and high cost, the present disclosure broadens the frequency band width of the amplitude-preserved pre-stack shot gather data and refines the depth-domain velocity model, fully utilizes all effective width frequency band information of the pre-stack data, and then reduces the bin size and increases the coverage times through pre-stack interpolation wavefield reconstruction, and applies the traveling wave separation to optimize the imaging condition in the reverse time migration imaging process to obtain a stacked seismic data volume, and applies the excision and diffusion filtering method on the extracted imaging gathers to double reduce the low wavenumber reverse time migration noise, comprehensively improve the ability to accurately depict the details of complex geological targets and the signal-to-noise ratio, realize wideband high-resolution reverse time migration processing, the applied frequency band is more than doubled compared with the conventional method, and finally efficiently and accurately improves the ability to accurately depict complex geological targets without increasing the computational complexity and cost.

[0129] The embodiments of the present disclosure have been described above. The above description is exemplary and not exhaustive, and is also not limited to the disclosed embodiments. Many modifications and variations are obvious to those of ordinary skill in the art without departing from the scope and spirit of the described embodiments. The choice of terms used herein is intended to best explain the principles of the embodiments, practical applications, or improvements to technologies in the market, or to enable other ordinary skilled persons in the art to understand the embodiments disclosed herein.

Claims

1. A prestack reverse time migration imaging method for broadband seismic waves, characterized in that, Including: Obtaining the original prestack shot gather data, depth-domain velocity model, and source data in the study area; Performing denoising, spherical spreading compensation, surface-consistent energy compensation, and deconvolution processing on the original prestack shot gather data to obtain a low-frequency-preserved, fidelity-preserved, wide-band prestack shot gather data; the deconvolution processing includes a combination of surface-consistent deconvolution and predictive deconvolution; the interpolation processing includes moveout correction processing, prestack interpolation processing, and reverse moveout correction processing; Using the depth-domain velocity model to optimize the parameters of the reverse time migration algorithm to obtain the reverse time migration algorithm with optimized parameters, the method including: setting the time difference accuracy of the reverse time migration numerical discretization algorithm to be greater than or equal to a first predetermined value, and the spatial difference accuracy to be greater than or equal to a second predetermined value; Under the condition of the difference accuracy, using the depth-domain velocity model to perform reverse time migration parameter optimization processing to obtain the reverse time migration algorithm with optimized parameters, the method including: obtaining the minimum velocity in the depth-domain velocity model; determining the number of grid nodes occupied by each seismic wave wavelength under the difference accuracy; setting the maximum bin and minimum depth step of the migration; according to the number of grid nodes, the minimum velocity, the maximum bin, and the minimum depth step, determining the maximum migration frequency fmax1 allowed for reverse time migration in the spatial direction and the maximum migration frequency fmax2 allowed for reverse time migration in the depth direction; according to the maximum migration frequency allowed for reverse time migration in the spatial direction and the maximum migration frequency allowed for reverse time migration in the depth direction, determining the maximum frequency parameter used for wide-band reverse time migration; where the maximum frequency parameter fmax = 2 * fmax1, and when the calculation result fmax >= fmax2, then fmax = fmax2; Before using the depth-domain velocity model to optimize the parameters of the reverse time migration algorithm, it further includes: performing prestack interpolation processing on the low-frequency-preserved, fidelity-preserved, wide-band prestack shot gather data in the common shot domain, and performing anti-spherical spreading compensation processing on the prestack shot gather data after the interpolation processing; Using the low-frequency-preserved, fidelity-preserved, wide-band prestack shot gather data, depth-domain velocity model, and source data as input data, and performing reverse time migration using the reverse time migration algorithm with optimized parameters to obtain a reverse time migration data volume; Extracting an imaging gather from the reverse time migration data volume; Performing excision stacking and low wavenumber noise removal processing on the imaging gather to obtain a reverse time migration imaging result.

2. The prestack reverse time migration imaging method for wideband seismic waves according to claim 1, wherein Before obtaining the depth-domain velocity model, it further includes: establishing a depth-domain velocity model; The method for establishing the depth-domain velocity model includes: Obtaining an isotropic velocity model vel or an anisotropic parameter model; Performing grid tomography global inversion on the isotropic velocity model vel or anisotropic parameter model to obtain a depth-domain velocity model; Wherein, the depth step of the grid tomography global inversion is within a first threshold, and the inversion grid step in the X and Y directions is less than a first predetermined value and greater than 2 bin sizes.

3. The wide-band seismic wave prestack reverse time migration imaging method according to claim 1, characterized in that: The first predetermined value is 2nd order, and the second predetermined value is 8th order.

4. The broadband seismic wave pre-stack reverse time migration imaging method according to claim 3, characterized in that: The maximum migration frequency fmax1 allowed for the spatial direction reverse time migration = velmin / (2.5 * binmax); where velmin is the minimum velocity in the depth domain velocity model, and binmax is the maximum bin for migration; The maximum migration frequency fmax2 allowed for the depth direction reverse time migration = velmin / (2.5 * depthmin), where depthmin is the minimum depth step.

5. The broadband seismic wave pre-stack reverse time migration imaging method according to claim 1, characterized in that The method of using the low-frequency-preserved and fidelity-preserved broadband pre-stack shot gather data, the depth domain velocity model, and the source data as input data, and performing reverse time migration using the optimized parameter reverse time migration algorithm to obtain the reverse time migration data volume includes: Using the optimized parameter reverse time migration algorithm, in the depth domain velocity model, performing wavefield continuation on the low-frequency-preserved and fidelity-preserved broadband pre-stack shot gather data and the source data to generate the receiver wavefield of the pre-stack shot gather data and the source wavefield, and performing up-going wave and down-going wave wavefield separation processing on the receiver wavefield of the pre-stack shot gather data and the source wavefield at each moment; According to the separated wavefields, determining the reverse time imaging condition at each time step for each imaging position; Accumulating the reverse time imaging conditions at all time steps corresponding to each imaging position to obtain the imaging point value at this imaging position, and all the imaging point values constitute the reverse time migration data volume; The method of performing wavefield continuation on the source wave data to generate the source wavefield includes: According to the maximum frequency parameter fmax and the band-limited wavelet function formula, obtaining a band-limited wavelet as the source; Using the discretized seismic wave equation in the reverse time migration algorithm, in the depth domain velocity model, generating the source wavefield through forward wavefield continuation; where the band-limited wavelet function formula includes: ； In the formula, is the upper limit frequency of the band-limited wavelet, that is, fmax, is the lower limit frequency of the band-limited wavelet, and is not equal to 0, and t is time.

6. The broadband seismic wave pre-stack reverse time migration imaging method according to claim 5, characterized in that, The method of determining the reverse time imaging condition at each time step for each imaging position according to the separated wavefields includes: According to the separated wavefields, using the reverse time imaging condition determination formula to determine the reverse time imaging condition at each time step for each position in the imaging space; Among them, the reverse time imaging condition determination formula includes: I = R 上 *S 上 +R 下 *S 下 ; Where: R 上 is the up-going wave of the detector point wavefield, S 上 is the up-going wave of the source wavefield, R 下 is the down-going wave of the detector point wavefield, S 下 is the down-going wave of the source wavefield.

7. The broadband seismic wave pre-stack reverse time migration imaging method according to claim 1, characterized in that The method of extracting the imaging gather in the reverse time migration data volume includes: In the reverse time migration data volume, extracting the data information of each shot trace header to form a common image point super gather; According to the set number of azimuths, starting offset, and offset interval, reconstructing the data traces of the imaging point super gather to obtain the imaging gather.

8. The broadband seismic wave pre-stack reverse time migration imaging method according to claim 1, wherein The method of mute stacking and low wavenumber noise removal processing includes: In the full coverage position of the study area, picking up a mute function on the imaging gather at any position, and applying the mute function to the imaging gather to obtain the imaging gather after removing the stretched and distorted wavefield; Accumulating the values of the imaging gather after removing the stretched and distorted wavefield along the numerical values in the same depth direction to obtain the reverse time migration stacked data volume; Using the diffusion filtering method to perform relative amplitude suppression processing on the low wavenumber noise of the reverse time migration stacked data volume to obtain the reverse time migration imaging result; The specific formula of the diffusion filtering method includes: ， ， ； Wherein: is the iterative filtering coefficient, N is the total number of iterative filtering times, t represents the current iteration number, and U 0 is the initial migration stack data volume, and U t+1 is the migration stack data volume after the (t + 1)-th iteration.

9. A prestack reverse time migration imaging device for broadband seismic waves, characterized in that, Including: An acquisition unit for acquiring the original prestack shot gather data, depth-domain velocity model, and source data of the study area; A prestack shot gather data processing unit for denoising, spherical divergence compensation, surface-consistent energy compensation, and deconvolution processing of the original prestack shot gather data to obtain a low-frequency-preserved, fidelity-preserved, wide-band prestack shot gather data; the deconvolution processing includes a combination of surface-consistent deconvolution and predictive deconvolution; the interpolation processing includes moveout correction processing, prestack interpolation processing, and reverse moveout correction processing; A parameter optimization unit for optimizing the reverse time migration algorithm parameters using the depth-domain velocity model to obtain an optimized reverse time migration algorithm, and the method includes: setting the time difference accuracy of the reverse time migration numerical discretization algorithm to be greater than or equal to a first predetermined value, and the spatial difference accuracy to be greater than or equal to a second predetermined value; under the condition of the difference accuracy, using the depth-domain velocity model to perform reverse time migration parameter optimization processing to obtain an optimized reverse time migration algorithm, and the method includes: obtaining the minimum velocity in the depth-domain velocity model; determining the number of grid nodes occupied by each seismic wave wavelength under the difference accuracy; setting the maximum bin and minimum depth step of the migration; according to the number of grid nodes, the minimum velocity, the maximum bin, and the minimum depth step, determining the maximum migration frequency fmax1 allowed for reverse time migration in the spatial direction and the maximum migration frequency fmax2 allowed for reverse time migration in the depth direction; according to the maximum migration frequency allowed for reverse time migration in the spatial direction and the maximum migration frequency allowed for reverse time migration in the depth direction, determining the maximum frequency parameter adopted for wide-band reverse time migration; where the maximum frequency parameter fmax = 2 * fmax1, and when the calculation result fmax >= fmax2, then fmax = fmax2; before using the depth-domain velocity model to optimize the reverse time migration algorithm parameters, it also includes: performing prestack interpolation processing on the low-frequency-preserved, fidelity-preserved, wide-band prestack shot gather data in the common shot point domain, and performing anti-spherical divergence compensation processing on the prestack shot gather data after the interpolation processing; A reverse time migration unit for using the low-frequency-preserved, fidelity-preserved, wide-band prestack shot gather data, depth-domain velocity model, and source data as input data, and performing reverse time migration using the optimized reverse time migration algorithm to obtain a reverse time migration data volume; An imaging gather extraction unit for extracting imaging gathers from the reverse time migration data volume; An imaging result generation unit for performing excision stacking and low wavenumber noise removal processing on the imaging gathers to obtain a reverse time migration imaging result.

Images