Viscoelastic prestack depth migration method and system based on path-dependent equivalent Q value

Abstract

The application provides a viscoelastic prestack depth migration method and system based on a path-dependent equivalent Q value, and particularly relates to the technical field of seismic exploration. The method comprises the following steps: collecting prestack seismic data of a target work area; performing anti-Q filtering on the prestack seismic data to obtain a constant Q scanning gather varying with a moveout group; performing constant Q scanning on the constant Q scanning gather according to a preset time window to determine a constant Q value corresponding to each time window; establishing a path-dependent equivalent Q field based on the constant Q values of all moveout groups; and performing viscoelastic prestack depth migration based on the path-dependent equivalent Q field to obtain an imaging result. The method does not depend on time-domain Q field modeling and travel time calculation, can greatly improve the efficiency and accuracy of depth-domain Q field modeling, and can significantly improve the imaging resolution and accuracy when the constructed depth-domain Q field is used for viscoelastic prestack depth migration imaging.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of seismic exploration, in particular to a viscoelastic prestack depth migration method and system based on path-dependent equivalent Q value. BACKGROUND

[0002] Viscoelastic prestack depth migration is the most effective method for high-resolution seismic imaging, and Q field modeling is a key step in viscoelastic prestack depth migration. Traditional Q value estimation relies on VSP logging, interwell data and other transmission wave information, or converts Q value from velocity field through empirical formula. However, due to insufficient data coverage and low signal-to-noise ratio, it is difficult to obtain large-scale three-dimensional high-precision non-uniform Q field by using these methods. It is feasible to obtain Q value by using surface seismic data, but the tuning phenomenon caused by thin layer stacking will seriously interfere with the main wave waveform, making it difficult to distinguish whether the change of the waveform is caused by intrinsic attenuation of the medium or thin layer tuning. If this interference is not eliminated, the Q value obtained by using the spectral amplitude may be incorrect, although some research has tried to solve this problem, but the thin layer tuning interference still needs to be effectively considered when estimating Q value.

[0003] Therefore, the viscoelastic prestack depth migration technology introduces the concept of equivalent Q value, and the local window constant Q scanning can effectively avoid the thin layer tuning problem, laying a foundation for establishing the depth domain Q field. The depth domain Q field modeling method in the prior art is to first establish a time domain equivalent Q value model, and then map it to the depth domain through time-depth conversion. When calculating the travel time table, the Q value of each point on the ray path can be obtained according to the existing real ray path, so as to perform subsequent viscoelastic prestack depth migration. It can be seen that this method needs to perform viscoelastic prestack time migration to obtain the time domain equivalent Q field before performing viscoelastic prestack depth migration, and the time domain Q value modeling and imaging ray time-depth conversion will inevitably introduce certain errors, especially in the area with complex underground structure and severe lateral velocity variation, which is easy to affect the accuracy of the depth domain Q field, thereby causing the resolution of the imaging result to be low. In addition, the viscoelastic prestack depth migration using the calculated virtual travel time technology also faces many challenges such as travel time table calculation and storage. SUMMARY

[0004] In view of the above deficiencies of the prior art, the purpose of the present application is to provide a viscoelastic prestack depth migration method and system based on path-dependent equivalent Q value, which aims to solve the calculation and storage problems caused by the need for time domain Q field modeling and virtual travel time calculation in the existing viscoelastic prestack depth migration method.

[0005] In order to achieve the above purpose, the first aspect of the present application provides a viscoelastic prestack depth migration method based on path-dependent equivalent Q value, mainly comprising the following steps:

[0006] Collecting prestack seismic data of the target work area;

[0007] performing inverse Q filtering on the pre-stack seismic data to obtain depth domain constant Q scan gathers varying with offset groups;

[0008] performing constant Q scanning on the constant Q scan gathers according to preset time windows to determine constant Q values of the offset groups corresponding to each of the time windows;

[0009] establishing a path-dependent equivalent Q field based on the constant Q values of all the offset groups;

[0010] performing viscoelastic pre-stack depth migration based on the path-dependent equivalent Q field to obtain an imaging result.

[0011] Optionally, the performing constant Q scanning on the constant Q scan gathers according to preset time windows to determine constant Q values of the offset groups corresponding to each of the time windows comprises:

[0012] dividing the constant Q scan gathers according to preset time windows to determine a plurality of Q values corresponding to each of the time windows;

[0013] calculating stacking spectra corresponding to all the Q values in each window, and determining constant Q values of the offset groups corresponding to each of the time windows based on the stacking spectra.

[0014] Optionally, the establishing a path-dependent equivalent Q field based on the constant Q values of all the offset groups comprises:

[0015] establishing a depth domain Q value model based on a parabolic relationship between Q values and offsets;

[0016] determining a best parameter combination of the depth domain Q value model corresponding to each time window based on the constant Q values of the offset groups;

[0017] establishing a path-dependent equivalent Q field based on the best parameter combination.

[0018] Optionally, the determining a best parameter combination of the depth domain Q value model corresponding to each time window based on the constant Q values of the offset groups comprises:

[0019] constructing an objective function based on constant Q spectra and path-dependent Q spectra corresponding to each time window;

[0020] solving an optimal solution of the objective function corresponding to each time window by using a preset optimization algorithm based on the constant Q values of the offset groups;

[0021] determining a best parameter combination of the depth domain Q value model corresponding to each time window based on the optimal solution.

[0022] Optionally, the target function is constructed based on the normal Q spectrum and the path-dependent Q spectrum corresponding to each time window, comprising:

[0023] The center frequency, the main low-frequency component, the frequency width and the spectrum area of the normal Q spectrum and the path-dependent Q spectrum corresponding to each time window are calculated.

[0024] Based on the center frequency, the main low-frequency component, the frequency width and the spectrum area, the superposition spectrum index and the partial offset spectrum index corresponding to the normal Q spectrum and the path-dependent Q spectrum are constructed respectively.

[0025] Based on all the superposition spectrum indexes and the partial offset spectrum indexes, the target function is constructed.

[0026] Optionally, the optimal solution of the target function corresponding to each time window is solved by using a preset optimization algorithm based on the normal Q value of the offset group, comprising:

[0027] Based on the normal Q value of the offset group, the parameters in the depth domain Q value model are initialized to obtain an initial value.

[0028] Based on the superposition spectrum and the partial offset spectrum corresponding to the normal Q spectrum, the initial solution of the target function is obtained.

[0029] Based on the initial value and the initial solution, the optimal solution of the target function corresponding to each time window is solved by using a probability-based global optimization algorithm.

[0030] Optionally, the path-dependent equivalent Q field is established based on the best parameter combination, comprising:

[0031] Based on the best parameter combination corresponding to each time window, the depth domain Q value model is solved to obtain the best Q value corresponding to each time window.

[0032] All the best Q values are interpolated and smoothed to establish the path-dependent equivalent Q field.

[0033] The second aspect of the present application provides a viscoelastic prestack depth migration system based on path-dependent equivalent Q value, comprising:

[0034] An information acquisition module is configured to collect prestack seismic data of a target work area.

[0035] A normal Q scan gather generation module is configured to perform inverse Q filtering on the prestack seismic data to obtain depth domain normal Q scan gathers varying with offset groups.

[0036] A constant Q value generation module is configured to perform constant Q scanning on the constant Q scanning gathers according to a preset time window, and determine a constant Q value of each offset group corresponding to the time window;

[0037] A path-dependent equivalent Q field construction module is configured to establish a path-dependent equivalent Q field based on the constant Q values of all the offset groups;

[0038] A viscoelastic prestack depth migration module is configured to perform viscoelastic prestack depth migration based on the path-dependent equivalent Q field, and obtain an imaging result.

[0039] The third aspect of the present application provides an intelligent terminal, which comprises a memory, a processor, and a viscoelastic prestack depth migration program based on a path-dependent equivalent Q value stored in the memory and executable on the processor, and the viscoelastic prestack depth migration program based on the path-dependent equivalent Q value implements the steps of any one of the above viscoelastic prestack depth migration methods based on the path-dependent equivalent Q value when executed by the processor.

[0040] The fourth aspect of the present application provides a computer readable storage medium, which stores a viscoelastic prestack depth migration program based on a path-dependent equivalent Q value, and the viscoelastic prestack depth migration program based on the path-dependent equivalent Q value implements the steps of any one of the above viscoelastic prestack depth migration methods based on the path-dependent equivalent Q value when executed by a processor.

[0041] Compared with the prior art, the present application has the following advantages:

[0042] The viscoelastic prestack depth migration method based on a path-dependent equivalent Q value provided by the present application is based on the parabolic correlation between Q value and offset distance, and by introducing parameters such as main frequency, main low frequency component and frequency width to describe the spectral characteristics, the offset distance is divided into multiple offset distance groups, the problem of determining the optimal parabolic parameter combination is transformed into a multi-parameter optimization problem by comparing the viscoelastic prestack depth migration result spectrum of the constant Q value of the offset distance group, and the fast simulated annealing method is used for solving, so as to directly establish a depth domain Q field model. This method does not depend on time domain Q field modeling, reduces the error and a large amount of manual and computing resource cost caused by time domain Q value modeling, and does not depend on travel time calculation, which can greatly improve the efficiency and accuracy of depth domain Q field modeling. The depth domain Q field constructed by the method can be used for viscoelastic prestack depth migration imaging, and can significantly improve the imaging resolution and accuracy. BRIEF DESCRIPTION OF DRAWINGS

[0043] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the drawings needed to be used in the embodiments or prior art description will be briefly introduced. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained by those skilled in the art without any creative effort.

[0044] Figure 1 Flow chart of the viscoelastic pre-stack depth migration method based on path-dependent equivalent Q value of the present application;

[0045] Figure 2 Contrastive schematic diagram of the relationship between different depth migration distances of the present application;

[0046] Figure 3 Schematic diagram of low-frequency, main-frequency and high-frequency parameters of the present application;

[0047] Figure 4 Schematic diagram of the conventional pre-stack depth migration result of the 415th line in the imaging space of the work area of the present application;

[0048] Figure 5 Path-dependent depth-domain Q field contour map of the present application;

[0049] Figure 6 Schematic diagram of the viscoelastic pre-stack depth migration result obtained from the path-dependent equivalent Q field of the present application;

[0050] Figure 7 Schematic diagram of the viscoelastic pre-stack depth migration result of the 415th line in the imaging space of the work area (1500m-4000m) of the present application;

[0051] Figure 8 Schematic diagram of the viscoelastic pre-stack depth migration result of the 415th line in the imaging space of the work area (2500m-5000m) of the present application;

[0052] Figure 9 Schematic diagram of the viscoelastic pre-stack depth migration result of the 415th line in the imaging space of the work area (4200m-6000m) of the present application;

[0053] Figure 10 Schematic diagram of the viscoelastic pre-stack depth migration result of the 415th line in the imaging space of the work area (4200m-6000m) of the present application; Figure 7-Figure 9 Contrastive diagram of the pre-stack depth migration result spectrum corresponding to the window of the present application;

[0054] Figure 11 Schematic diagram of the viscoelastic pre-stack depth migration system module based on path-dependent equivalent Q value of the present application;

[0055] Figure 12 Schematic diagram of the intelligent terminal structure of the present application. DETAILED DESCRIPTION

[0056] In the following description, for purposes of explanation and not limitation, specific details are set forth such as particular architectures, techniques, etc. in order to provide a thorough understanding of the embodiments of the present application. However, it will be apparent to those skilled in the art that the present application can be practiced in other embodiments that depart from these specific details. In other instances, detailed descriptions of well-known methods, devices, circuits, and

[0057] It is to be understood that the terminology "includes", "has", "holds", "contains" and / or "comprising", "comprised of", "comprising", "comprises" when used in this specification and in the following claims, specifies the presence of stated features, integers, steps, operations, elements, and / or components but do not preclude the presence or addition of one or more other features, integers, steps, operations, elements, components, and / or groups thereof.

[0058] It is also to be understood that the terminology used in the description of the present application is for the purpose of describing particular embodiments only and is not intended to be limiting of the present application. As used in this description and the appended claims, the singular forms "a", "an", and "the" are intended to include the plural forms as well, unless the context clearly indicates otherwise.

[0059] It will be further understood that the terms "comprises" and / or "comprising", when used in this specification and in the following claims, specify the presence of stated features, integers, steps, operations, elements, and / or components, but do not preclude the presence or addition of one or more other features, integers, steps, operations, elements, components, and / or groups thereof.

[0060] The technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings of the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, rather than all the embodiments of the present application. Based on the embodiments of the present application, all other embodiments obtained by a person of ordinary skill in the art without creative effort fall within the scope of protection of the present application.

[0061] Many specific details are set forth in the following description in order to provide a thorough understanding of the present application. However, persons of ordinary skill in the art will realize that the present application can be practiced without many of these details. In other instances, well-known structures and functions have not been described in detail in order to avoid obscuring the present application. The description is made for the purpose of illustration only and is not intended to be limiting.

[0062] In order to solve the problem that the existing depth domain Q value modeling needs to rely on time domain Q field modeling, so as to improve the viscoelastic prestack depth migration Q field modeling precision, the application provides a viscoelastic prestack depth migration method based on path correlation equivalent Q value, mainly based on the characteristics that Q value and offset distance are parabolic correlation, main frequency, main low frequency component and frequency width and other parameters for describing spectral characteristics are introduced, offset distance is divided into multiple offset distance groups, by comparing with the viscoelastic prestack depth migration result spectrum of the constant Q value of the offset distance group, the problem of determining the best parabolic parameter combination is converted into a multi-parameter optimization problem, and the fast simulated annealing method is used for solving, so as to directly establish the depth domain Q value model. The method directly establishes the depth domain Q field model without relying on the time domain Q field and imaging ray time-depth conversion for modeling, can greatly improve the efficiency and precision of the depth domain Q field modeling, and is beneficial to improving the resolution of the viscoelastic prestack depth migration imaging.

[0063] The embodiment of the application provides a viscoelastic prestack depth migration method based on path correlation equivalent Q value, which is deployed on a computer, a server or other electronic equipment, applied to scenes related to geological structure such as seismic exploration, stratum exploration, oil and gas exploration, and is aimed at the case of viscoelastic prestack depth migration based on depth domain Q field analysis. Specifically, as shown in the figure, Figure 1 The steps of the method include:

[0064] Step S100: collecting prestack seismic data of a target work area;

[0065] Specifically, a suitable exploration area is selected as the target work area according to the exploration purpose (such as oil and gas exploration, geological structure research, etc.). A reasonable seismic line is designed according to the exploration target, including line direction, interval, length, etc. Seismic waves are excited on the line at a certain interval and in a certain way, and high-precision seismic recording instruments are used to record the prestack seismic data received by the geophone. These data contain rich underground medium information.

[0066] Further, the prestack seismic data can be pretreated, such as denoising, filtering, dynamic correction and static correction, before subsequent step processing, so as to improve the quality of the prestack seismic data and reduce the error in subsequent processing.

[0067] Step S200: performing inverse Q filtering on the prestack seismic data to obtain depth domain constant Q scan gathers varying with offset distance groups;

[0068] Specifically, inverse Q filtering is a seismic data processing method for compensating the Q attenuation effect in seismic records, making the spectrum of seismic records flat, thereby improving the resolution and interpretation ability of seismic images. In the processing of pre-stack seismic data, the application of inverse Q filtering can improve the data quality. The quality factor (Q value) is an important parameter for measuring the degree of attenuation of seismic waves in underground media, which is used to reflect the amplitude change caused by the propagation of rays from one point to another. By accumulating the Q value attenuation effect on the entire propagation path, the Q value effect reflecting the amplitude and phase change from the shot point to the imaging point can be obtained. By analyzing the attenuation characteristics of pre-stack seismic data, such as the change of amplitude with propagation distance, frequency spectrum attenuation, etc., the Q value of the underground medium can be estimated. According to the estimated Q value range, a series of Q values are selected as the input parameters for subsequent inverse Q filtering, and the selected Q values can cover the possible attenuation range of the underground medium, so as to comprehensively evaluate the influence of different Q values on seismic data.

[0069] When the reflection interface is horizontal, different traces in different common shot gathers on the survey line can always be found, which are all from a common point on the underground interface, which is called common depth point or common reflection point. The corresponding record of the common depth reflection point is called CDP (or CRP) gather. The CRP gather is the result of rearranging seismic data according to the reflection point position, which can more clearly show the characteristics of the underground reflection interface.

[0070] In seismic exploration, the Q value (quality factor) is closely related to the attenuation of seismic waves, and the attenuation is affected by the characteristics of underground media (such as lithology, porosity, fluid content, etc.). Ray tracing is performed using an actual depth domain velocity model, and ray tracing is performed at depths of 4000m, 5000m and 7000m respectively at the center of the earth's surface. The relationship between offset h and Q value Q h at different depths is compared as shown in Figure 2 , wherein Figure 2 a), b) and c) in represent the simulation results at depths of 4000m, 5000m and 7000m respectively. It can be observed that the offset represents the distance from the source to the receiver. Therefore, based on the parabolic relationship between the Q value and the offset, the offset is divided into multiple groups of partial offsets, and the Q value corresponding to different groups of partial offsets is different, and it is related to the parabolic relationship between the partial offset of the group and the partial offset of the group. For the inverse Q filtering result corresponding to each Q value, the corresponding common reflection point (CRP) gather can be generated, and all CRP gathers corresponding to different Q values are stacked to obtain a depth domain constant Q scanning gather varying with offset group.

[0071] Step S300: performing a constant-Q scan on the constant-Q scan gather according to a preset time window, and determining a constant-Q value of an offset group corresponding to each time window;

[0072] Specifically, the constant-Q scanning gather is divided according to preset time windows, and several Q values ​​corresponding to each time window are determined; the superimposed spectrum corresponding to all the Q values ​​in each window is calculated, and based on the superimposed spectrum, the constant Q value of the offset group corresponding to each time window is determined to obtain the Q value changes in different time windows, which helps to understand the changes in underground lithology and its corresponding physical properties.

[0073] Based on the characteristics of the subsurface geology in the target area, a series of time windows are preset. The selection of time window size should take into account the complexity of the subsurface geology and the resolution of the seismic data. Excessively large windows may obscure variations in geological details, while too small windows may lead to inaccurate Q-value estimates due to insufficient data. These windows can be equally spaced or unevenly spaced based on geological characteristics such as geological horizons. The constant-Q scan gather data is divided into preset time windows, with each time window containing seismic trace data within a set time range. For each time window, all seismic trace data are extracted. For each seismic trace, the corresponding Q-value is calculated using a Q-value estimation algorithm such as the spectral ratio method or waveform fitting, based on the characteristics of amplitude attenuation and phase distortion during seismic wave propagation. Then, for each time window, the spectra corresponding to the Q-values ​​calculated for all seismic traces within the corresponding offset group in that window are superimposed to generate a superimposed spectrum, which reflects the average effect of the seismic wave attenuation characteristics within that time window. Finally, based on the superimposed spectrum, a constant Q-value that best describes the seismic wave attenuation characteristics within that time window is fitted or estimated to represent the average attenuation characteristics of the subsurface medium within that time window.

[0074] In this embodiment, constant-Q scanning gathers are effectively divided and processed to obtain the attenuation characteristics of underground media in different time windows, providing an important basis for seismic exploration interpretation and underground geological structure analysis.

[0075] Step S400: establishing a path-dependent equivalent Q field based on the constant Q values ​​of all the offset groups;

[0076] Specifically, based on Figure 2 The offset h and Q value Q at different depths are shown h Based on the parabolic relationship between Q value and offset, this embodiment divides the offset into multiple component offsets. Different component offsets correspond to different Q values, and they are parabolically correlated with the offsets of the group. For a single time window, a depth domain Q value model is constructed, namely:

[0077] Qh = a h 2 + b h + Q0 (1)

[0078] where h represents a set of offset distances, a, b and Q0 are parameters of the parabolic model.

[0079] Then, based on the offset distances h and Q values at different depths obtained by ray tracing using the actual depth domain velocity model, the best parameter combination (a, b, Q0) of the depth domain Q value model corresponding to each time window is determined. best Considering that the attenuation characteristics of seismic waves in the underground are related to not only the depth but also the specific propagation path (for example, there is additional attenuation when passing through different lithological interfaces), it is necessary to analyze the influence of the path on the Q value. Based on the best parameter combination of each time window and the actual propagation path of the seismic wave, a path-dependent equivalent Q field is constructed, which not only considers the influence of the depth on the Q value but also considers the modulation effect of the path on the Q value, which is beneficial to provide more accurate geological information for seismic exploration. However, since the number of parameter combination (a, b, Q0) samples obtained by ray tracing using the actual depth domain velocity model is small, it is not easy to find the best parameter combination of the depth domain Q value model.

[0080] Since for each selected time window, the imaging result of the path-dependent Q field is affected by the stacking spectrum and the spectrum corresponding to each set of offset distances. Therefore, in order to improve the resolution of the viscoelastic prestack depth migration imaging result corresponding to the path-dependent Q field, an evaluation index M (Metric) is constructed based on the correlation characteristics of the stacking spectrum and the spectrum corresponding to each set of offset distances to improve the characteristics of the stacking spectrum and the spectrum corresponding to each set of offset distances, thereby compensating the viscoelastic prestack depth migration imaging result. Therefore, based on the constant Q spectrum and the path-dependent Q spectrum corresponding to each time window, a target function is constructed to convert the best parameter combination problem of the depth domain Q value model into a multi-parameter optimization problem with the target function M. Since the values of the target function M corresponding to different parameter combinations are different. Then, the global optimal solution of the target function M corresponding to each time window can be solved by using a preset optimization algorithm; based on the global optimal solution, the best parameter combination of the depth domain Q value model corresponding to each time window is determined, thereby establishing the path-dependent equivalent Q field.

[0081] In a preferred embodiment, the construction of the target function based on the constant Q spectrum and the path-dependent Q spectrum corresponding to each time window comprises:

[0082] The center frequency, the main low frequency component, the frequency width and the spectrum area of the normal Q spectrum and the path-related Q spectrum corresponding to each time window are calculated; based on the center frequency, the main low frequency component, the frequency width and the spectrum area, the superimposed spectrum index and the split offset spectrum index corresponding to the normal Q spectrum and the path-related Q spectrum respectively are constructed; and based on all the superimposed spectrum indexes and the split offset spectrum indexes, the objective function is constructed.

[0083] Specifically, the spectrum characteristics are described by four characteristics of the central frequency Cf (Central Frequency), the main low frequency component f min (FrequencyMinimum), the frequency width f bw (Frequency Bandwidth) and the spectrum area SA (Spectrum Area), and the evaluation indexes are constructed by using the spectrum characteristics.

[0084] The expression of the central frequency CF (Central Frequency) is as follows:

[0085]

[0086] Wherein, F is the frequency value, and S is the amplitude value corresponding to the frequency.

[0087] Since the viscoelastic pre-stack depth migration imaging result is compensated, the main frequency is improved, the high frequency signal is recovered, and the low frequency signal is also lost, but the ideal situation is to keep the low frequency basically unchanged and improve the high frequency, which means that the high and low frequency parts of the spectrum are improved as much as possible, and intuitively, it is to ensure that the flatness TF (TopFlatness) of the spectrum is as large as possible. The expression of the flatness TF of the spectrum is as follows:

[0088]

[0089] Wherein, N represents the number of frequency sampling points.

[0090] From the above formula, it can be seen that the flatness TF of the spectrum can only reflect the relative flatness, and cannot reflect the absolute frequency width. Directly using TF to measure the spectrum characteristics is not rigorous. In fact, the frequency width f bw (Frequency Bandwidth) represents the difference between the main high frequency component f max (Frequency Maximum) and the main low frequency component f min (Frequency Minimum). And CF and TF can be used to determine f min , f max , f min and f maxRespectively on both sides of CF, ideally both high and low frequencies within the bandwidth are enhanced, which means f min to f max are as large as possible. Assuming a given threshold threshold, extending to both sides of the center frequency CF, when TF reaches the threshold, the left and right boundaries of the spectrum are f min and f max , and f bw =f max -f min +1. Figure 3 The low frequency, main frequency and high frequency parameter diagram is shown.

[0091] The expression of the spectrum area SA is as follows:

[0092] SA=∑S (4)

[0093] When both low and high frequencies are enhanced, the spectrum area SA is also enhanced.

[0094] For each selected time window, the imaging result of the path-dependent Q field corresponds to the spectrum, which not only enhances the relevant characteristics of the stacked spectrum, but also needs to enhance the spectrum corresponding to each offset distance group. The common Q spectrum indicators are: stacked spectrum characteristics S CF 0 , and S SA 0 , and the expression of the characteristics of the offset distance spectrum: CF 0 , f min 0 , f bw 0 and SA 0 , wherein S in the stacked spectrum characteristics is the abbreviation of stack; The indicators of the path-dependent Q are: stacked spectrum characteristics S CF * , and S SA * , and the expression of the characteristics of the offset distance spectrum: CF * , f min * , f bw * and SA * . The number of offset distances is n, and the above characteristics of the offset distance spectrum are the average values of the spectrum characteristics corresponding to n different offset distances, that is:

[0095]

[0096] Finally, based on the above constant Q spectrum and path-dependent Q spectrum corresponding to each time window, the expression of the evaluation index M (Metric) is as follows:

[0097]

[0098] It can be seen that the evaluation index M is composed of the above eight parts, each part is the ratio of the path-dependent Q index and the constant Q index, and the value of f min is the smaller the better, and the values of the remaining items are the larger the better, so the correlation term of f min is set in inverse proportion, and the remaining terms are set in direct proportion. For each term, the value is less than 1, indicating that the corresponding feature is worse than the constant Q index, and vice versa, the larger the value, the better the corresponding feature than the constant Q index.

[0099] In a preferred embodiment, the constant Q value based on the offset distance group is used to solve the optimal solution of the objective function corresponding to each time window by using a preset optimization algorithm, including:

[0100] Based on the constant Q value of the offset distance group, the parameters in the depth domain Q value model are initialized to obtain an initial value; based on the constant Q spectrum corresponding to the superposition spectrum and the split offset distance spectrum, the initial solution of the objective function is obtained; based on the initial value and the initial solution, the optimal solution of the objective function corresponding to each time window is solved by using a probability-based global optimization algorithm.

[0101] Specifically, since the simulated annealing algorithm is a probability-based global optimization algorithm, it finds the global optimal solution by simulating the actual metal annealing process in physics, and is suitable for solving complex optimization problems. Therefore, the simulated annealing algorithm is used for optimization in this embodiment, and the algorithm process is as follows:

[0102] 1) Parameter initialization. Set the initial temperature T, the temperature threshold of each iteration dT, and the initial parameters P0 and the initial solution M, and the initial temperature is related to the size of the solution space, T needs to be high enough to accept more poor solutions in the early stage and seek more possible solutions;

[0103] 2) Loop until the stop condition is met. A small random perturbation is added to the neighborhood of the current parameter P to obtain a new parameter P * , and a new solution M * corresponding to the new parameter is calculated * , the difference ΔM=M * -M between the current parameter P and the new parameter P * is calculated, and when ΔM>0, that is, the new solution is better than the current solution, the new parameter P -ΔM / T is accepted, and when ΔM<0, that is, the new solution is worse than the current solution, the new parameter P *For the suboptimal solution is not necessarily abandoned, but a certain degree of acceptance to jump out of the local optimal solution.

[0104] 3) According to T i = T i-1 × dT temperature update, constantly updated to the value of the temperature T is less than the threshold Thres then stop iteration to get the global optimal solution.

[0105] Through the above steps, the simulated annealing algorithm can effectively explore the solution space and avoid falling into a local optimal solution, gradually approaching the global optimal solution, and realizing global optimization.

[0106] In this embodiment, the indicators used to construct the target function M can be divided into stack spectrum indicators and offset distance spectrum indicators. The stack spectrum indicators directly reflect the resolution and signal-to-noise ratio of the final imaging result, and the offset distance spectrum indicators, as a reference for the final imaging result, should have a relatively lower weight than the stack spectrum indicators. In actual application, the importance of the stack spectrum indicators and the offset distance spectrum indicators can be adjusted by setting the weights of the four stack spectrum indicators and the weights of the four offset distance spectrum indicators, so that the resolution and signal-to-noise ratio of the viscoelastic prestack depth migration imaging result reach the expectation. For example, the weights of each stack spectrum indicator can be set to be the same (such as 1.2 / 8), and the weights of each offset distance spectrum indicator can be set to be the same (such as 0.8 / 8), or the same or different weights can be set for each or some of them according to needs.

[0107] The constant Q full offset stack spectrum is obtained from the CRP gather corresponding to the constant Q value of all offset distance groups, and the constant Q stack spectrum feature S CF 0 、 and S SA 0 ; According to the time window, the offset distance is divided into n groups. In order to ensure that multi-parameter optimization can be performed, n should be greater than or equal to 3, and considering that the data in each group of offset distance should not be too small, the value of n should not be too large, and n is usually preferably 5, 6 or 7. Assuming that each group of offset distance contains m Q values, that is, the offset distance is divided into n groups from the CRP gather corresponding to the m Q values, a total of n x m groups of offset distance data, the offset distance stack spectrum is calculated, and n x m groups of offset distance spectrum features CF * , f min * , f bw * and SA * , and according to f max * , the corresponding high cut-off frequency F3 * is determined. The constant Q offset distance spectrum feature in the n x m group offset distance spectrum feature is: CF 0 , fmin 0 , f bw 0 and SA 0 .

[0108] The initial value P0=(a, b, Q0) of parameter optimization is determined, and the Q value corresponding to each offset distance group can be obtained from the initial parameters: (h1, Q1), (h2, Q2)…(h n , Q n ), the constant Q stack result is obtained by full offset stack of the CRP gather corresponding to the constant Q, and the path-dependent Q stack data is obtained by first determining the Q value corresponding to the offset distance, then taking the offset distance data from the CRP gather of the Q value, and then obtaining the path-dependent Q stack spectrum by stacking all the partial offset distance data, and then determining the stacking spectrum feature S CF * , and S SA * After obtaining all the parameters above, the initial solution M0 can be determined. Based on the above simulated annealing algorithm process, after the initial parameters P0 and the initial solution M are determined, the initial temperature T, the temperature reduction coefficient dT of each iteration, and the temperature threshold T are set to globally optimize the target function M, obtain the global optimal solution, and thus determine the best parameter combination (a, b, Q0) corresponding to the time window.

[0109] Similarly, the optimal parameter combination corresponding to each time window is solved according to the above process, and the parabolic parameter combination of the window Q value and offset distance is obtained, so as to obtain a set of optimal Q values corresponding to each time window. Then, the optimal Q values are interpolated and smoothed to establish the path-dependent equivalent Q field. Specifically, the optimal Q values corresponding to all time windows are sorted according to their corresponding depths or positions to form an initial Q value data point set. According to the characteristics of the Q value data and the geological conditions of the study area, a suitable interpolation method is selected, such as linear interpolation, nearest neighbor interpolation, bilinear interpolation, bicubic interpolation, Kriging interpolation, etc. Among them, Kriging interpolation is a statistical interpolation method, which can consider the spatial correlation and variability of data points, and is particularly suitable for geological data interpolation. The selected interpolation method is used to interpolate the Q value data point set to fill in the missing Q values between the data points and form a continuous Q value surface. Then, a smoothing method such as moving average, Gaussian filtering, median filtering, etc. is selected to smooth the Q value surface formed after interpolation, and the Q value data is filtered or weighted averaged, etc. to eliminate noise and irregularities, reduce noise and irregularities in the Q value data, and make the Q field model smoother. Finally, the path-dependent equivalent Q field is established to accurately represent the Q value changes of the underground medium at different depths and paths, thereby accurately reflecting the attenuation characteristics of the seismic wave when passing through different media or complex geological structures.

[0110] Further, the optimal frequency bandwidth and the optimal main low-frequency component can also be obtained from the global optimal solution bw = f max -f min +1, and the optimal main high-frequency component , i.e. the high cut-off frequency F3, corresponding to different sub-offset distances is obtained. The high cut-off frequency under different sub-offset distances is mapped to the spatial or depth domain as input data. Using geological statistics, machine learning or physical models, a mathematical model of the F3 field is constructed according to the spatial distribution of the high cut-off frequency to reflect the characteristics of the F3 field changing with offset distance and underground medium.

[0111] Step S500: performing viscoelastic prestack depth migration based on the path-dependent equivalent Q field to obtain an imaging result.

[0112] Specifically, using the path-dependent equivalent Q field for viscoelastic prestack depth migration processing can well restore the detailed information of the underground geological structure, compensate for the energy attenuation of the seismic wave when propagating along different paths, and help to restore the attenuated high-frequency component, thereby improving the resolution and accuracy of imaging.

[0113] Further, the imaging results of viscoelastic prestack depth migration can be analyzed to identify underground structure features such as faults, lithologic interfaces, etc. Based on the imaging results and geological interpretation, reservoir prediction and evaluation are performed to evaluate the distribution range, thickness, physical property and other parameters of the reservoir, thereby providing reliable basis for subsequent oil and gas exploration and development.

[0114] In summary, compared with the existing viscoelastic prestack depth migration technology relying on time-domain Q field modeling and travel time calculation, in actual imaging process, due to the large travel time table, it is difficult to directly store in the memory, and the Q field calculated by the travel time also faces this problem, therefore, in actual application, all the travel times and Q value effects are not completely calculated, but a multiple optimization interpolation method is used to reduce the size to ensure that it can be stored in the GPU memory for imaging. Obviously, converting the time-domain Q field to the depth domain and multiple optimization interpolation will cause errors, thus the resolution and accuracy of the final migration imaging are low, and the interpretability of the inversion result is poor. The method of directly performing depth-domain Q field modeling designed in the present application is based on the parabolic correlation between the Q value and the offset, by introducing parameters such as main frequency, main low-frequency component and frequency width to describe the spectral characteristics, by comparing with the spectral results of the viscoelastic prestack depth migration with constant Q value, the problem of determining the optimal parabolic parameter combination is converted into a multi-parameter optimization problem, and a fast simulated annealing method is used to solve it, thereby directly establishing the depth-domain Q field model. This method does not rely on time-domain Q field modeling, reduces the errors caused by time-domain Q value modeling and the cost of a large amount of manual and computing resources, and does not rely on travel time calculation, which can greatly improve the efficiency and accuracy of depth-domain Q field modeling. The depth-domain Q field constructed for viscoelastic prestack depth migration imaging can significantly improve the imaging resolution and accuracy, and is conducive to improving the industrial application prospect of viscoelastic prestack depth migration.

[0115] The effectiveness of the method of the present application is verified by the description of the examples in combination with the accompanying drawings.

[0116] The example data is three-dimensional actual data of a work area in eastern China, the imaging space has 1131 sampling points along the line direction with a sampling interval of 10 meters, 450 sampling points along the vertical line direction with a sampling interval of 20 meters, and 800 sampling points in the depth direction with a sampling interval of 10 meters. The imaging range is defined as 0.7 km to 11.3 km with an imaging interval of 10 meters, and the output profile depth is 8 km with a depth sampling interval of 3 meters.

[0117] Firstly, the prestack seismic data is migrated by using the conventional prestack depth migration technology to obtain the conventional migration result, such as Figure 4The figure shows the conventional pre-stack depth migration result of the 415th line of the imaging space of the work area. Then, the inverse Q filtering and data preparation, the constant Q scanning, the determination of all time window parameter combinations by optimizing the objective function using the simulated annealing algorithm, and the interpolation smoothing of the parameter combinations are sequentially performed according to the steps of the method of the application, to obtain the path-dependent Q field, as shown in Figure 5 The figure shows the path-dependent depth domain Q field contour map, the horizontal coordinate represents the CDP gather, the vertical coordinate represents the depth, and the contour inside the figure represents the offset, Figure 5 The figures a) to f) respectively show the Q fields corresponding to different offset groups, the offsets are 50m-783m, 783m-1516m, 1516m-2249m, 2249m-2982m, 2982m-3715m and 3715m-4450m, respectively, and the figures show the path-dependent Q field distribution map, and it can be seen that the Q fields of different offsets are non-uniformly distributed. The viscoelastic pre-stack depth migration is performed using the path-dependent Q field, as shown in Figure 6 The figure shows the viscoelastic pre-stack depth migration result obtained based on the path-dependent equivalent Q field, and it can be seen that the result has higher resolution than Figure 4 The figure shows the conventional pre-stack depth migration result. Further, the viscoelastic pre-stack depth migration of different depths is respectively performed using the conventional pre-stack depth migration method and the viscoelastic pre-stack depth migration method of the application, Figure 7 The figure shows the viscoelastic pre-stack depth migration (1500m-4000m) result of the 415th line of the imaging space of the work area, Figure 7 The figures a) and b) respectively show the imaging result of the conventional method and the imaging result of the method of the application; Figure 8 The figure shows the viscoelastic pre-stack depth migration (2500m-5000m) result of the 415th line of the imaging space of the work area, Figure 8 The figures a) and b) respectively show the imaging result of the conventional method and the imaging result of the method of the application; Figure 9 The figure shows the viscoelastic pre-stack depth migration (4200m-6000m) result of the 415th line of the imaging space of the work area, Figure 9 The figures a) and b) respectively show the imaging result of the conventional method and the imaging result of the method of the application. From the imaging results from the shallow layer to the deep layer, Figure 7-Figure 9 The figure shows the imaging results from the shallow layer to the deep layer, and compared with the conventional method, the resolution of the imaging result obtained by the method of the application is obviously improved, and the events are more continuous and well separated. Finally, Figure 10 The figures a) to c) respectively show Figure 7-Figure 9 The figure shows the spectrum comparison of the pre-stack depth migration result of the corresponding window, and it can be seen that the effective frequency width of the viscoelastic migration result of the method of the application is significantly improved compared with the conventional method. The above simulation experiment fully proves the effectiveness and feasibility of the method of the application.

[0118] The figure shows the imaging result of the 415th line of the imaging space of the work area, Figure 11As shown, corresponding to the viscoelastic pre-stack depth migration method based on path-dependent equivalent Q value, the embodiment of the present application also provides a viscoelastic pre-stack depth migration system based on path-dependent equivalent Q value, the viscoelastic pre-stack depth migration system based on path-dependent equivalent Q value comprises:

[0119] The information acquisition module 111 is configured to collect pre-stack seismic data of a target work area.

[0120] The constant Q scanning trace set generation module 112 is configured to perform inverse Q filtering on the pre-stack seismic data to obtain depth domain constant Q scanning trace sets varying with offset distance groups.

[0121] The constant Q value generation module 113 is configured to perform constant Q scanning on the constant Q scanning trace sets according to a preset time window to determine constant Q values of offset distance groups corresponding to each time window.

[0122] The path-dependent equivalent Q field construction module 114 is configured to establish a path-dependent equivalent Q field based on the constant Q values of all the offset distance groups.

[0123] The viscoelastic pre-stack depth migration module 115 is configured to perform viscoelastic pre-stack depth migration based on the path-dependent equivalent Q field to obtain an imaging result.

[0124] Specifically, in the embodiment, the specific functions of the viscoelastic pre-stack depth migration system based on path-dependent equivalent Q value can also be referred to the corresponding description in the viscoelastic pre-stack depth migration method based on path-dependent equivalent Q value, which will not be described here.

[0125] Based on the above embodiment, the present application further provides an intelligent terminal, and a principle block diagram thereof can be as shown. Figure 12 The intelligent terminal comprises a processor, a memory, a network interface and a display screen connected through a system bus. The processor of the intelligent terminal is configured to provide calculation and control capabilities. The memory of the intelligent terminal comprises a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system and a viscoelastic pre-stack depth migration program based on path-dependent equivalent Q value. The internal memory provides an environment for the operating system and the viscoelastic pre-stack depth migration program based on path-dependent equivalent Q value in the non-volatile storage medium. The network interface of the intelligent terminal is configured to communicate with external terminals through network connection. The viscoelastic pre-stack depth migration program based on path-dependent equivalent Q value, when executed by the processor, implements the steps of any one of the viscoelastic pre-stack depth migration methods based on path-dependent equivalent Q value. The display screen of the intelligent terminal can be a liquid crystal display screen or an electronic ink display screen.

[0126] Those skilled in the art can understand that, Figure 12The principle block diagram shown in the figure is only a block diagram of part of the structure related to the scheme of the present application, and does not constitute a limitation on the intelligent terminal to which the scheme of the present application is applied. The specific intelligent terminal can include more or fewer components than those shown in the figure, or combine certain components, or have a different component arrangement.

[0127] In one embodiment, an intelligent terminal is provided, which comprises a memory, a processor, and a path-dependent equivalent Q-value-based viscoelastic pre-stack depth migration program stored in the memory and executable on the processor. The path-dependent equivalent Q-value-based viscoelastic pre-stack depth migration program, when executed by the processor, implements the steps of any path-dependent equivalent Q-value-based viscoelastic pre-stack depth migration method provided in the embodiments of the present application.

[0128] The embodiments of the present application also provide a computer readable storage medium, which stores a path-dependent equivalent Q-value-based viscoelastic pre-stack depth migration program. The path-dependent equivalent Q-value-based viscoelastic pre-stack depth migration program, when executed by a processor, implements the steps of any path-dependent equivalent Q-value-based viscoelastic pre-stack depth migration method provided in the embodiments of the present application.

[0129] It should be understood that the sequence numbers of the steps in the above embodiments do not mean the order of execution. The execution order of the processes should be determined according to their functions and inherent logic, and should not constitute any limitation on the implementation process of the embodiments of the present application.

[0130] Those skilled in the art can clearly understand that, for the convenience and brevity of description, only the division of the above functional units and modules is taken as an example. In actual application, the above functions can be completed by different functional units and modules according to needs, that is, the internal structure of the device is divided into different functional units or modules to complete all or part of the above-described functions. The functional units and modules in the embodiments can be integrated in one processing unit, or each unit can exist physically, or two or more units can be integrated in one unit. The integrated unit can be realized in the form of hardware or in the form of a software functional unit. In addition, the specific names of the functional units and modules are only for easy distinction, and do not limit the protection scope of the present application. The specific working process of the units and modules in the system can refer to the corresponding process in the foregoing method embodiments, which will not be repeated here.

[0131] In the above embodiments, the description of each embodiment has its own emphasis. The parts not described or recorded in detail in a certain embodiment can be referred to the relevant description of other embodiments.

[0132] Those skilled in the art can understand that the units and algorithm steps of each example described in combination with the embodiments disclosed herein can be realized in electronic hardware or a combination of computer software and electronic hardware. Whether the functions are performed in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different ways to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of the present application.

[0133] In the embodiments provided by the present application, it should be understood that the disclosed apparatus / terminal device and method can be implemented in other ways. For example, the apparatus / terminal device embodiments described above are merely schematic, and the division of the above modules or units is merely a logical function division, and there can be another division in actual implementation. For example, a plurality of units or components can be combined or integrated into another system, or some features can be ignored or not executed.

[0134] The above-described embodiments are merely used to illustrate the technical solutions of the present application, rather than limit them; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand; the technical solutions recorded in the foregoing embodiments can be modified, or some technical features can be replaced by equivalents; and these modifications or replacements do not deviate from the spirit and scope of the corresponding technical solutions, and should be included in the protection scope of the present application.

Claims

1. A viscoelastic prestack depth migration method based on path-dependent equivalent Q values, characterized by, The method comprises the following steps: Collecting pre-stack seismic data of a target work area; Performing inverse Q filtering on the pre-stack seismic data to obtain depth domain constant Q scan gathers varying with offset distance groups; Performing constant Q scanning on the constant Q scan gathers according to preset time windows to determine constant Q values of offset distance groups corresponding to each time window; Based on the constant Q values of all the offset distance groups, establishing a path-dependent equivalent Q field; Based on the path-dependent equivalent Q field, performing viscoelastic pre-stack depth migration to obtain imaging results; The step of establishing a path-dependent equivalent Q field based on the constant Q values of all the offset distance groups comprises: Based on a parabolic relationship between Q values and offset distances, establishing a depth domain Q value model; Based on the constant Q values of the offset distance groups, determining a best parameter combination of the depth domain Q value model corresponding to each time window; Based on the best parameter combination, establishing a path-dependent equivalent Q field; The step of determining a best parameter combination of the depth domain Q value model corresponding to each time window based on the constant Q values of the offset distance groups comprises: Based on constant Q spectra and path-dependent Q spectra corresponding to each time window, constructing an objective function; Based on the constant Q values of the offset distance groups, solving an optimal solution of the objective function corresponding to each time window by using a preset optimization algorithm; Based on the optimal solution, determining a best parameter combination of the depth domain Q value model corresponding to each time window; The step of constructing an objective function based on constant Q spectra and path-dependent Q spectra corresponding to each time window comprises: Calculating center frequencies, main low-frequency components, frequency widths and spectral areas of the constant Q spectra and the path-dependent Q spectra corresponding to each time window; Based on the center frequencies, the main low-frequency components, the frequency widths and the spectral areas, constructing superposition spectral indexes and split offset distance spectral indexes corresponding to the constant Q spectra and the path-dependent Q spectra respectively; Based on all the superposition spectral indexes and the split offset distance spectral indexes, constructing an objective function.

2. The method of claim 1, wherein the path-dependent equivalent Q-value based viscoelastic prestack depth migration method is characterized by, The step of performing constant Q scanning on the constant Q scan gathers according to preset time windows to determine constant Q values of offset distance groups corresponding to each time window comprises: Dividing the constant Q scan gathers according to preset time windows to determine a plurality of Q values corresponding to each time window; Calculating superposition spectra of all the Q values in each window, and based on the superposition spectra, determining constant Q values of offset distance groups corresponding to each time window.

3. The method of claim 1, wherein the path-dependent equivalent Q-value based viscoelastic prestack depth migration method is characterized by, The step of solving an optimal solution of the objective function corresponding to each time window by using a preset optimization algorithm based on the constant Q values of the offset distance groups comprises: Initializing parameters in the depth domain Q value model to obtain initial values; Based on superposition spectra and split offset distance spectra corresponding to the constant Q spectra, obtaining an initial solution of the objective function; Based on the initial values and the initial solution, solving an optimal solution of the objective function corresponding to each time window by using a probability-based global optimization algorithm.

4. The method of claim 1, wherein, The step of establishing a path-dependent equivalent Q field based on the best parameter combination comprises: Based on the best parameter combination corresponding to each time window, solving the depth domain Q value model to obtain a best Q value corresponding to each time window; Interpolating and smoothing all the optimal Q values to build a path-dependent equivalent Q field.

5. A viscoelastic prestack depth migration system based on path-dependent equivalent Q values, characterized by, The system comprises: An information acquisition module configured to collect prestack seismic data of a target work area; A constant Q scan gather generation module configured to perform inverse Q filtering on the prestack seismic data to obtain constant Q scan gathers in depth domain varying with offset groups; A constant Q value generation module configured to perform constant Q scanning on the constant Q scan gathers according to preset time windows to determine constant Q values of offset groups corresponding to each of the time windows; A path-dependent equivalent Q field construction module configured to build a path-dependent equivalent Q field based on the constant Q values of all the offset groups; A viscoelastic prestack depth migration module configured to perform viscoelastic prestack depth migration based on the path-dependent equivalent Q field to obtain an imaging result; The building of the path-dependent equivalent Q field based on the constant Q values of all the offset groups comprises: building a depth domain Q value model based on a parabolic relationship between Q values and offsets; determining, based on the constant Q values of the offset groups, an optimal parameter combination of the depth domain Q value model corresponding to each of the time windows; building the path-dependent equivalent Q field based on the optimal parameter combination; The determination of the optimal parameter combination of the depth domain Q value model corresponding to each of the time windows based on the constant Q values of the offset groups comprises: constructing an objective function based on constant Q spectra and path-dependent Q spectra corresponding to each of the time windows; solving, based on the constant Q values of the offset groups, an optimal solution of the objective function corresponding to each of the time windows by using a preset optimization algorithm; determining, based on the optimal solution, the optimal parameter combination of the depth domain Q value model corresponding to each of the time windows; The construction of the objective function based on the constant Q spectra and the path-dependent Q spectra corresponding to each of the time windows comprises: calculating center frequencies, main low-frequency components, frequency widths and spectral areas of the constant Q spectra and the path-dependent Q spectra corresponding to each of the time windows; constructing, based on the center frequencies, the main low-frequency components, the frequency widths and the spectral areas, superposition spectral indexes and split-offset spectral indexes corresponding to the constant Q spectra and the path-dependent Q spectra respectively; constructing the objective function based on all the superposition spectral indexes and the split-offset spectral indexes.

6. A smart terminal, characterized by The intelligent terminal comprises a memory, a processor, and a viscoelastic prestack depth migration program based on path-dependent equivalent Q values stored on the memory and executable on the processor, and the viscoelastic prestack depth migration program based on path-dependent equivalent Q values, when executed by the processor, implements the steps of the viscoelastic prestack depth migration method based on path-dependent equivalent Q values as claimed in any one of claims 1-4.

7. A computer-readable storage medium, characterized in that, The computer readable storage medium stores a viscoelastic prestack depth migration program based on path-dependent equivalent Q values, and the viscoelastic prestack depth migration program based on path-dependent equivalent Q values, when executed by the processor, implements the steps of the viscoelastic prestack depth migration method based on path-dependent equivalent Q values as claimed in any one of claims 1-4.

Images