Viscoelastic pre-stack depth migration method and system based on path-related equivalent Q value

Abstract

The invention provides a viscoelastic pre-stack depth migration method and system based on a path-dependent equivalent Q value, and particularly relates to the technical field of seismic exploration. The scheme 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 a constant Q scanning gather changing along with the offset group; performing constant Q scanning on the constant Q scanning gather according to preset time windows, and determining a constant Q value corresponding to each time window; establishing a path-related equivalent Q field based on the constant Q values of all the offset groups; and performing viscoelastic pre-stack depth migration based on the path-dependent equivalent Q field to obtain an imaging result. The scheme does not depend on time domain Q field modeling and travel time calculation, the efficiency and precision of depth domain Q field modeling can be greatly improved, the constructed depth domain Q field is used for viscoelastic pre-stack depth migration imaging, and the imaging resolution and accuracy can be remarkably improved.

Description

Technical Field

[0001] The present invention relates to the field of seismic exploration technology, and in particular to a viscoelastic prestack depth migration method and system based on path-related equivalent Q value. Background Art

[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 transmission wave information such as VSP logging and interwell data, or obtains Q value from velocity field through empirical formula conversion. However, due to factors such as insufficient data coverage and low signal-to-noise ratio, these methods are difficult to obtain large-scale three-dimensional high-precision non-uniform Q field. It is feasible to use surface seismic data to obtain Q value. However, the thin layer tuning phenomenon caused by thin layer stacking will seriously interfere with the main wave waveform, making it difficult to distinguish whether the waveform change is due to intrinsic attenuation of the medium or thin layer tuning. If this interference is not eliminated, the use of spectrum amplitude may result in an erroneous Q value. Although some studies have tried to solve this problem, thin layer tuning interference still needs to be effectively considered when estimating Q value.

[0003] Based on this, the viscoelastic prestack depth migration technology introduces the concept of equivalent Q value. The thin layer tuning problem can be effectively avoided by scanning with constant Q in a local window, laying the foundation for establishing the Q field in the depth domain. The depth domain Q field modeling method in the prior art is to first establish an equivalent Q value model in the time domain, and then map it to the depth domain through time-to-depth conversion. When obtaining 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. Time domain Q value modeling and imaging ray time-to-depth conversion will inevitably introduce certain errors, especially in areas with complex underground structures and drastic lateral velocity changes, which can easily affect the accuracy of the depth domain Q field, resulting in a low resolution of the imaging results. In addition, the prior art of using virtual travel time calculation technology to perform viscoelastic prestack depth migration will also face many challenges such as travel time table calculation and storage. Summary of the invention

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

[0005] In order to achieve the above object, the present invention provides a viscoelastic prestack depth migration method based on path-dependent equivalent Q value in a first aspect, which mainly includes the following steps:

[0006] Collect pre-stack seismic data in the target work area;

[0007] Performing inverse Q filtering on the pre-stack seismic data to obtain a depth domain constant Q scanning gather that varies with the offset group;

[0008] Performing a constant Q scan on the constant Q scan gathers according to a preset time window to determine a constant Q value of an offset group corresponding to each time window;

[0009] establishing a path-dependent equivalent Q field based on a constant Q value for all offset groups;

[0010] Viscoelastic prestack depth migration is performed based on the path-dependent equivalent Q field to obtain imaging results.

[0011] Optionally, performing a constant Q scan on the constant Q scan gathers according to a preset time window to determine a constant Q value of an offset group corresponding to each time window includes:

[0012] Dividing the constant-Q scan gathers according to preset time windows, and determining a number of Q values ​​corresponding to each of the time windows;

[0013] The superimposed frequency spectra corresponding to all the Q values ​​in each window are calculated, and based on the superimposed frequency spectra, the constant Q value of the offset group corresponding to each time window is determined.

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

[0015] Based on the parabolic relationship between Q value and offset, a depth domain Q value model is established;

[0016] Determining an optimal parameter combination of the depth domain Q value model corresponding to each time window based on the constant Q value of the offset group;

[0017] Based on the optimal parameter combination, a path-dependent equivalent Q field is established.

[0018] Optionally, determining the optimal parameter combination of the depth domain Q value model corresponding to each time window based on the constant Q value of the offset group includes:

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

[0020] Based on the constant Q value of the offset group, a preset optimization algorithm is used to solve the optimal solution of the objective function corresponding to each time window;

[0021] Based on the optimal solution, the best parameter combination of the depth-domain Q-value model corresponding to each time window is determined.

[0022] Optionally, constructing the objective function based on the constant Q spectrum and the path-related Q spectrum corresponding to each time window includes:

[0023] Calculate the center frequency, main low-frequency component, bandwidth and spectrum area of ​​the constant Q spectrum and path-related Q spectrum corresponding to each time window;

[0024] Based on the center frequency, the main low-frequency component, the bandwidth and the spectrum area, constructing superimposed spectrum indicators and offset spectrum indicators corresponding to the constant Q spectrum and the path-related Q spectrum respectively;

[0025] Based on all the superimposed spectrum indicators and the offset spectrum indicators, an objective function is constructed.

[0026] Optionally, the method of solving the optimal solution of the objective function corresponding to each time window by using a preset optimization algorithm based on the constant Q value of the offset group includes:

[0027] Initializing parameters in the depth-domain Q-value model based on the constant Q value of the offset group to obtain initial values;

[0028] Obtaining an initial solution of the objective function based on a superimposed spectrum and an offset spectrum corresponding to the constant-Q spectrum;

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

[0030] Optionally, establishing a path-dependent equivalent Q field based on the optimal parameter combination includes:

[0031] Solving the depth domain Q value model based on the optimal parameter combination corresponding to each of the time windows to obtain the optimal Q value corresponding to each of the time windows;

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

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

[0034] Information acquisition module, used to collect pre-stack seismic data of the target work area;

[0035] A constant Q scanning gather generation module is used to perform inverse Q filtering on the pre-stack seismic data to obtain a depth domain constant Q scanning gather that varies with the offset group;

[0036] A constant Q value generating module, used for 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;

[0037] A path-dependent equivalent Q field building module is used to build a path-dependent equivalent Q field based on the constant Q value of all the offset groups;

[0038] The viscoelastic prestack depth migration module is used to perform viscoelastic prestack depth migration based on the path-related equivalent Q field to obtain imaging results.

[0039] A third aspect of the present invention provides an intelligent terminal, which includes a memory, a processor, and a viscoelastic prestack depth migration program based on path-correlated equivalent Q values ​​stored in the memory and executable on the processor. When the viscoelastic prestack depth migration program based on path-correlated equivalent Q values ​​is executed by the processor, any step of the above-mentioned viscoelastic prestack depth migration method based on path-correlated equivalent Q values ​​is implemented.

[0040] A fourth aspect of the present invention provides a computer-readable storage medium, on which is stored a viscoelastic prestack depth migration program based on a path-correlated equivalent Q value, and when the viscoelastic prestack depth migration program based on a path-correlated equivalent Q value is executed by a processor, any step of the above-mentioned viscoelastic prestack depth migration method based on a path-correlated equivalent Q value is implemented.

[0041] Compared with the prior art, the beneficial effects of this solution are as follows:

[0042] The viscoelastic prestack depth migration method based on path-related equivalent Q value proposed in the present invention is based on the parabolic correlation between Q value and offset distance. By introducing parameters describing the spectrum characteristics such as main frequency, main low-frequency component and bandwidth, the offset distance is divided into multiple offset distance groups. By comparing with the spectrum of viscoelastic prestack depth migration results of constant Q value of the offset distance group, the problem of determining the best parabolic parameter combination is transformed into a multi-parameter optimization problem, and the fast simulated annealing method is used to solve it, so as to directly establish a depth domain Q field model. The method does not rely on time domain Q field modeling, reduces the error caused by time domain Q value modeling and a large amount of labor and computing resource costs, and does not rely on travel time calculation, which can greatly improve the efficiency and accuracy of depth domain Q field modeling. The constructed depth domain Q field is used for viscoelastic prestack depth migration imaging, which can significantly improve imaging resolution and accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] In order to more clearly illustrate the technical solutions in the embodiments of the present invention, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0044] Figure 1 It is a flow chart of the viscoelastic prestack depth migration method based on path-dependent equivalent Q value of the present invention;

[0045] Figure 2 It is a schematic diagram comparing the relationship between the offset distance and at different depths of the present invention;

[0046] Figure 3 It is a schematic diagram of low frequency, main frequency and high frequency parameters of the present invention;

[0047] Figure 4 This is a schematic diagram of conventional prestack depth migration results for line 415 of the work area imaging space of the present invention;

[0048] Figure 5 It is the path-related depth domain Q field contour map of the present invention;

[0049] Figure 6 A schematic diagram of the viscoelastic prestack depth migration result obtained from the path-dependent equivalent Q field of the present invention;

[0050] Figure 7 It is a schematic diagram of the results of viscoelastic prestack depth migration (1500 meters-4000 meters) of the 415th line of the imaging space of the work area of ​​the present invention;

[0051] Figure 8 It is a schematic diagram of the viscoelastic prestack depth migration result (2500 meters-5000 meters) of the 415th line of the imaging space of the work area of ​​the present invention;

[0052] Fig. 9 It is a schematic diagram of the results of the viscoelastic prestack depth migration (4200 meters - 6000 meters) of the 415th line of the imaging space of the work area of ​​the present invention;

[0053] Fig.10 For the present invention Figure 7-Figure 9 Spectral comparison of prestack depth migration results for the corresponding window;

[0054] Fig.11 It is a schematic diagram of a module of a viscoelastic prestack depth migration system based on path-dependent equivalent Q value of the present invention;

[0055] Fig.12 It is a schematic diagram of the structure of the intelligent terminal of the present invention. DETAILED DESCRIPTION

[0056] In the following description, specific details such as specific system structures, technologies, etc. are provided for the purpose of illustration rather than limitation, so as to provide a thorough understanding of the embodiments of the present invention. However, it should be clear to those skilled in the art that the present invention may be implemented in other embodiments without these specific details. In other cases, detailed descriptions of well-known systems, devices, circuits, and methods are omitted to prevent unnecessary details from obstructing the description of the present invention.

[0057] It should be understood that when used in this specification and the appended claims, the term "comprising" indicates the presence of described features, integers, steps, operations, elements and / or components, but does not exclude the presence or addition of one or more other features, integers, steps, operations, elements, components and / or groups thereof.

[0058] It should also be understood that the terms used in the present specification are only for the purpose of describing specific embodiments and are not intended to limit the present invention. As used in the present specification and the appended claims, unless the context clearly indicates otherwise, the singular forms "a", "an" and "the" are intended to include plural forms.

[0059] It should be further understood that the term "and / or" used in the present description and the appended claims refers to any and all possible combinations of one or more of the associated listed items, and includes these combinations.

[0060] The following is a clear and complete description of the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0061] In the following description, many specific details are set forth to facilitate a full understanding of the present invention, but the present invention may also be implemented in other ways different from those described herein, and those skilled in the art may make similar generalizations without violating the connotation of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0062] In order to solve the problem that the existing depth domain Q value modeling needs to rely on the time domain Q field modeling, so as to improve the accuracy of Q field modeling of viscoelastic prestack depth migration, the present invention proposes a viscoelastic prestack depth migration method based on path-related equivalent Q value, which is mainly based on the characteristic that the Q value is parabolically correlated with the offset distance, introduces parameters such as the main frequency, main low-frequency component and bandwidth to describe the spectrum characteristics, divides the offset distance into multiple offset distance groups, and compares the spectrum of the viscoelastic prestack depth migration results of the constant Q value of the offset distance group, transforms the problem of determining the best parabolic parameter combination into a multi-parameter optimization problem, and adopts the fast simulated annealing method to solve it, so as to directly establish the depth domain Q value model. The method of the present invention realizes the direct establishment of the depth domain Q field model without relying on the time domain Q field and the imaging ray time-to-depth conversion for modeling, which can greatly improve the efficiency and accuracy of the depth domain Q field modeling, and is conducive to improving the resolution of viscoelastic prestack depth migration imaging.

[0063] The embodiment of the present invention provides a viscoelastic prestack depth migration method based on path-related equivalent Q value, which is deployed on electronic devices such as computers and servers and applied to scenes related to geological structures such as seismic exploration, stratigraphic exploration, oil and gas exploration, etc., and is aimed at the situation of viscoelastic prestack depth migration based on depth domain Q field analysis. Specifically, Figure 1 As shown, the steps of the method in this embodiment include:

[0064] Step S100: collecting pre-stack seismic data of the target work area;

[0065] Specifically, according to the exploration purpose (such as oil and gas exploration, geological structure research, etc.), a suitable exploration area is selected as the target work area. According to the exploration goal, a reasonable seismic survey line is designed, including the survey line direction, spacing, length, etc. Seismic waves are excited on the survey line at a certain interval and in a certain manner, and the pre-stack seismic data received by the detector is recorded using a high-precision seismic recording instrument. These data contain rich underground medium information.

[0066] Furthermore, before the subsequent processing steps, the pre-stack seismic data may be pre-processed by denoising, filtering, dynamic correction and static correction to improve the quality of the pre-stack seismic data and reduce errors in subsequent processing.

[0067] Step S200: performing inverse Q filtering on the pre-stack seismic data to obtain a depth domain constant Q scanning gather that varies with the offset group;

[0068] Specifically, inverse Q filtering is a seismic data processing method used to compensate for 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 data quality. The quality factor (Q value) is an important parameter to measure the attenuation degree of seismic waves in the underground medium. It is used to reflect the amplitude change caused by the propagation of rays from one point to another. By accumulating all the Q value attenuation effects on the propagation path, the Q value effect reflecting the amplitude phase change from the shot detection 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, the attenuation of the frequency spectrum, 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 of the 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 impact of different Q values ​​on seismic data.

[0069] When the reflection interface is horizontal, different tracks can always be found in different common shot point gathers on the survey line. They all come from a common point on the underground interface, which is called a common depth point or common reflection point. The corresponding recording tracks with common depth reflection points form a common depth point (or common reflection point) gather called a CDP gather (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, which in turn is affected by the characteristics of the underground medium (such as lithology, porosity, fluid content, etc.). Ray tracing is performed using the actual depth domain velocity model. Ray tracing is performed at depths of 4,000 meters, 5,000 meters, and 7,000 meters at the center of the surface, respectively, to obtain the offset distance h and Q value Q at different depths. h The relationship comparison diagram is as follows: Figure 2 As shown. Among them, Figure 2 a), b) and c) in the figure represent the simulation results corresponding to the depths of 4000 meters, 5000 meters and 7000 meters respectively. It can be clearly observed that they are parabolically correlated, where the offset represents the distance from the source to the receiver. Therefore, based on the parabolic relationship between the Q value and the offset, this embodiment divides the offset into multiple components of offsets. Different component offsets correspond to different Q values, and they are parabolically correlated with the offsets of the group. For ease of understanding, this embodiment refers to each component of offset as an offset group. For each inverse Q filtering result corresponding to each Q value, a corresponding common reflection point (CRP) gather can be generated, and all CRP gathers corresponding to different Q values ​​are superimposed to obtain a depth domain constant Q scanning gather that varies with the 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 to determine a number of Q values ​​corresponding to each time window; 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 change of Q value in different time windows, which is helpful to understand the change of underground lithology and its corresponding physical properties.

[0073] According to the characteristics of the underground geological structure of the target work area, a series of time windows are preset. The selection of the time window size should take into account the complexity of the underground geological structure and the resolution of the seismic data. A window that is too large may obscure the changes in geological details, while a window that is too small may lead to inaccurate Q value estimation due to insufficient data. These windows can be divided at equal intervals or at non-equal intervals according to geological features such as geological layers. The constant Q scanning gather data is divided according to the preset time windows, and each time window contains seismic trace data within the set time range. For each time window, all seismic trace data are extracted. For each seismic trace, based on the characteristics of amplitude attenuation and phase distortion of seismic waves during propagation, the corresponding Q value is calculated using Q value estimation algorithms such as spectrum ratio method and waveform fitting method. Then, for each time window, the frequency spectra corresponding to the Q values ​​calculated in the offset group corresponding to all seismic traces in the window are superimposed to obtain a superimposed spectrum to reflect the average effect of the seismic wave attenuation characteristics in the time window. Finally, a constant Q value that can best describe the seismic wave attenuation characteristics in the time window is fitted or estimated based on the superimposed spectrum to represent the average attenuation characteristics of the underground medium in the time window.

[0074] In this embodiment, the constant Q scanning gathers are effectively divided and processed to obtain the attenuation characteristics of the underground medium 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, that is:

[0077] Qh =αh 2 +βh+Q0 (1)

[0078] Where h represents a component offset, α, β and Q0 are the parameters of the parabolic model.

[0079] Then, based on the offset h and Q values ​​at different depths obtained by ray tracing using the actual depth-domain velocity model, the optimal parameter combination (α, β, Q0) of the depth-domain Q value model corresponding to each time window is determined. best . Considering that the attenuation characteristics of seismic waves when propagating underground are not only related to the depth, but also to the specific propagation path (for example, there will be additional attenuation when passing through the interface of different lithologies), it is necessary to analyze the impact of the path on the Q value. Based on the optimal parameter combination of each time window and the actual propagation path of the seismic wave, a path-dependent equivalent Q field is constructed. This equivalent Q field not only considers the impact of depth on the Q value, but also considers the modulation effect of the path on the Q value, which is conducive to providing more accurate geological information for seismic exploration. However, due to the small number of sample parameter combinations (α, β, Q0) obtained by ray tracing using the actual depth domain velocity model, it is not easy to find the optimal parameter combination for the depth domain Q value model.

[0080] Since for each selected time window, the imaging result of the path-dependent Q field will be affected by the stacked spectrum and the spectrum corresponding to each group of offsets. Therefore, in order to improve the resolution of the viscoelastic prestack depth migration imaging result corresponding to the path-dependent Q field, this embodiment constructs an evaluation index M (Metric) based on the correlation characteristics of the stacked spectrum and the spectrum corresponding to each group of offsets to improve the characteristics of the stacked spectrum and the spectrum corresponding to each group of offsets, thereby compensating the viscoelastic prestack depth migration imaging result. Therefore, this embodiment constructs an objective function based on the constant Q spectrum and the path-dependent Q spectrum corresponding to each time window to convert the optimal parameter combination problem of the depth domain Q value model into a multi-parameter optimization problem with the objective function M. Since the values ​​of the objective function M corresponding to different parameter combinations are different. Then, the preset optimization algorithm can be used to solve the global optimal solution of the objective function M corresponding to each time window; based on the global optimal solution, the optimal parameter combination of the depth domain Q value model corresponding to each time window is determined, thereby establishing a path-dependent equivalent Q field.

[0081] In a preferred embodiment, constructing the objective function based on the constant Q spectrum and the path-related Q spectrum corresponding to each time window includes:

[0082] The center frequency, main low-frequency component, bandwidth and spectrum area of ​​the constant 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 bandwidth and the spectrum area, the superimposed spectrum index and the offset spectrum index corresponding to the constant Q spectrum and the path-related Q spectrum are constructed; based on all the superimposed spectrum indexes and the offset spectrum indexes, an objective function is constructed.

[0083] Specifically, the center frequency Cf (Central Frequency), the main low-frequency component f min (FrequencyMinimum), bandwidth f bw Four characteristics, such as frequency bandwidth and spectrum area SA (spectrum area), describe the spectrum characteristics, and these spectrum characteristics are used to construct evaluation indicators.

[0084] Among them, the expression of center frequency CF (Central Frequency) is as follows:

[0085]

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

[0087] Since the compensation of viscoelastic prestack depth migration imaging results will increase the main frequency and restore the high-frequency signal, it will often cause some loss of low-frequency signal. However, the ideal situation is to keep the low frequency basically unchanged and increase the high frequency, which means that both the high and low frequency parts of the spectrum are increased as much as possible. 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] Where N represents the number of frequency sampling points.

[0090] From the above formula, we can see that the flatness of the spectrum TF can only reflect the relative flatness, but not the absolute bandwidth. It is not rigorous to directly use TF to measure the spectrum characteristics. In fact, the bandwidth f bw (Frequency Bandwidth) represents the main high-frequency component f max (Frequency Maximum) and the main low-frequency component f min (Frequency Minimum). And using CF and TF, we can determine f min With f max , f min With f maxOn both sides of CF, the ideal situation is that both high and low frequencies within the bandwidth are improved, which also means that f min to f max The flatness TF is as large as possible. Assuming a given threshold, the spectrum is extended to both sides with the center frequency CF as the center. When TF reaches the threshold, the left and right boundaries of the spectrum are f min With f max , and satisfies f bw =f max -f min +1. Figure 3 Shown is a schematic diagram of low frequency, main frequency and high frequency parameters.

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

[0092] SA=∑S (4)

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

[0094] For each selected time window, the corresponding spectrum of the imaging result of the path-related Q field not only improves the correlation characteristics of the stacked spectrum, but also needs to improve the corresponding spectrum of the offset group. Common Q spectrum indicators are: stacked spectrum characteristics S CF 0 , and S SA 0 , and the characteristic expression of the offset spectrum: CF 0 、f min 0 、f bw 0 and SA 0 , where S in the superposition spectrum feature is the abbreviation of stack; the indicators of path-related Q are: superposition spectrum feature S CF * , and S SA * , and the characteristic expression of the offset spectrum: CF * 、f min * 、f bw * and SA * The number of sub-offsets is n, and the above sub-offset spectrum features are the average values ​​of the spectrum features corresponding to n different offsets, that is:

[0095]

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

[0097]

[0098] It can be seen that the evaluation index M is composed of the above eight items, each of which is the ratio of the path-related Q index to the constant Q index, except f min The smaller the better, the values ​​of the other items are the larger the better, so f min The relevant items are set to be inversely proportional, and the remaining items are set to be directly proportional. For each item, its value is less than 1, indicating that the corresponding feature is worse than the constant Q index. Conversely, the larger its value, the better the corresponding feature is than the constant Q index.

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

[0100] Based on the constant Q value of the offset group, the parameters in the depth domain Q value model are initialized to obtain initial values; based on the superimposed spectrum and the offset spectrum corresponding to the constant Q 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 using a probability-based global optimization algorithm.

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

[0102] 1) Parameter initialization. Set the initial temperature T, the cooling coefficient dT and the temperature threshold for each iteration, as well as the initial parameter P0 and the initial solution M. The initial temperature is related to the size of the solution space. T needs to be high enough to ensure that more inferior solutions can be accepted in the early stage and more solutions can be sought;

[0103] 2) Loop until the stop condition is met. A random small perturbation is made in the neighborhood of the current parameter P to obtain a new parameter P * , and calculate the new solution M corresponding to the new parameters * , calculate the current parameter P and the new parameter P * The difference in objective function ΔM=M * -M, when ΔM>0, that is, the new solution is better than the current solution, accept the new parameter P * , when ΔM<0, that is, the new solution is inferior to the current solution, then with probability e -ΔM / T Accept new parameter P *, the inferior solution is not necessarily abandoned, but accepted to a certain extent in order to escape from the local optimal solution.

[0104] 3) According to T i =T i-1 ×dT to update the temperature, and continuously update it until the value temperature T is less than the threshold Thres, then stop the iteration and obtain the global optimal solution.

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

[0106] In this embodiment, the indexes used to construct the objective function M can be divided into stack spectrum indexes and offset spectrum indexes. The stack spectrum index directly reflects the resolution and signal-to-noise ratio of the final imaging result. The offset spectrum index as a reference for the final imaging result should have a relatively lower weight than the stack spectrum index. In practical applications, the importance of the stack spectrum index and the offset spectrum index can be adjusted by setting the weights of the four stack spectrum indexes and the weights of the four offset spectrum indexes, so that the resolution and signal-to-noise ratio of the viscoelastic prestack depth migration imaging result can reach the expected value. For example, the weights of the stack spectrum indexes are all set to be the same (for example, 1.2 / 8), and the weights of the offset spectrum indexes are all set to be the same (for example, 0.8 / 8). The same or different weights can also be set for each item or several of them as needed.

[0107] The constant Q full offset stack spectrum is obtained from the CRP gathers corresponding to the constant Q values ​​of all offset groups, and the constant Q stack spectrum feature S is determined. CF 0 , and S SA 0 ; Divide the offsets into n groups according to the time window. To ensure that multi-parameter optimization can be performed, n should be greater than or equal to 3. In addition, considering that the data in each offset group is not too small, the value of n should not be too large. It is usually preferred that n be 5, 6 or 7. Assume that each offset group contains m Q values, that is, the offsets are divided into n groups from the CRP gathers corresponding to the m Q values, with a total of n×m offset data. The offset stacking spectrum is calculated to obtain the n×m offset spectrum feature CF * 、f min * 、f bw * and SA * , and according to f max * Determine the corresponding high cutoff frequency F3 * The constant Q value in the n×m component offset spectrum characteristics corresponds to the constant Q offset spectrum characteristics: CF 0 、fmin 0 、f bw 0 and SA 0 .

[0108] Determine the initial value of parameter optimization P0 = (α, β, Q0), and the Q value corresponding to each offset group can be obtained from the initial parameters: (h1, Q1), (h2, Q2)… (h n , Q n ), the constant Q stacking result is obtained by stacking the CRP gathers corresponding to the constant Q at full offsets, and the path-dependent Q stacking data is obtained by first determining the Q value corresponding to the offset, then taking the offset data from the CRP gathers of the Q value, and then superimposing all the offset data to obtain the path-dependent Q stacking spectrum, and then determining the stacking spectrum characteristics S of the offsets. CF * , and S SA * , after obtaining all the above parameters, the initial solution M0 can be determined. Based on the above simulated annealing algorithm process, after determining the initial parameters P0 and the initial solution M, by setting the initial temperature T, the cooling coefficient dT of each iteration and the temperature threshold T, the objective function M is globally optimized to obtain the global optimal solution, thereby determining the optimal parameter combination (α, β, Q0) corresponding to the time window.

[0109] By analogy, the optimal parameter combination corresponding to each time window is solved according to the above process, and the depth domain Q value model is solved to obtain the parabolic parameter combination of the window Q value and the offset distance, so as to obtain a set of optimal Q values ​​corresponding to each time window. Then all the optimal Q values ​​are interpolated and smoothed to establish a path-related 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 set of Q value data points. 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 an interpolation method based on statistics. Because it can consider the spatial correlation and variability of data points, it is particularly suitable for the interpolation of geological data. Using the selected interpolation method, the Q value data point set is interpolated to fill the missing Q values ​​between the data points to form a continuous Q value surface. Then, smoothing methods such as moving average, Gaussian filtering, and median filtering are used to smooth the Q-value surface formed after interpolation, and the Q-value data are filtered or weighted averaged to eliminate noise and irregularities, reduce noise and irregularities in the Q-value data, and make the Q field model smoother. Finally, a path-related equivalent Q field is established to accurately characterize the Q-value changes of underground media at different depths and paths, thereby accurately reflecting the attenuation characteristics of seismic waves when passing through different media or through complex geological structures.

[0110] Furthermore, the optimal bandwidth can be obtained according to the global optimal solution. and optimal main low-frequency components According to the formula f bw =f max -f min +1, get the optimal main high-frequency components corresponding to different offsets That is, the high cutoff frequency F3. The high cutoff frequencies at different offsets are used as input data and mapped to the space or depth domain. Using methods such as geostatistics, machine learning or physical models, according to the spatial distribution law of the high cutoff frequency, a mathematical model of the F3 field is constructed to reflect the characteristics of the F3 field changing with the offset 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 to perform viscoelastic prestack depth migration processing can well restore detailed information on underground geological structures, compensate for the energy attenuation of seismic waves when propagating along different paths, help restore attenuated high-frequency components, and improve the resolution and accuracy of imaging.

[0113] Furthermore, the imaging results of viscoelastic prestack depth migration can be analyzed to identify underground structural features, such as faults and lithologic interfaces. Reservoir prediction and evaluation can be performed based on the imaging results and geological interpretation. The distribution range, thickness, physical properties and other parameters of the reservoir can be evaluated to provide a reliable basis for subsequent oil and gas exploration and development.

[0114] In summary, compared with the existing viscoelastic prestack depth migration technology that relies on time domain Q field modeling and travel time calculation, in the actual imaging process, the travel time table is too large to be directly stored in the memory. The Q field obtained by travel time calculation also faces this problem. Therefore, in practical applications, the travel time and Q value effects of all points are not completely calculated, but the multiple optimization interpolation method is used to reduce its 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, so the final migration imaging resolution and accuracy are low, making the inversion result less interpretable. The method designed by the present invention for directly modeling the depth domain Q field is based on the parabolic correlation between the Q value and the offset distance. By introducing parameters such as the main frequency, the main low-frequency component and the bandwidth that describe the spectrum characteristics, and by comparing with the spectrum of the viscoelastic prestack depth migration result of the constant Q value, 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 to solve it, so as to directly establish the depth domain Q field model. This method does not rely on time-domain Q-field modeling, which reduces the errors caused by time-domain Q-value modeling and a large amount of labor and computing resource costs. It also does not rely on travel time calculations, and can greatly improve the efficiency and accuracy of depth-domain Q-field modeling. The constructed depth-domain Q-field is used for viscoelastic prestack depth migration imaging, which can significantly improve imaging resolution and accuracy, and is conducive to improving the industrial application prospects of viscoelastic prestack depth migration.

[0115] The effectiveness of the method of the present invention is verified by describing examples in conjunction with the accompanying drawings.

[0116] The example data is the actual 3D data of a work area in eastern China. The data imaging space has 1131 sampling points along the survey line direction with a sampling interval of 10 meters, 450 sampling points along the vertical survey 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, the imaging interval is 10 meters, the output profile depth direction is 8 km, and the depth sampling interval is 3 meters.

[0117] First, the conventional prestack depth migration technology is used to migrate the prestack seismic data to obtain conventional migration results, such as Figure 4The figure shows the conventional prestack depth migration result of the 415th line of the imaging space of the work area. Then, according to the steps of the method of the present invention, the inverse Q filtering and data preparation, constant Q scanning, the use of simulated annealing algorithm to optimize the objective function to determine all time window parameter combinations, and the interpolation and smoothing of the parameter combinations are performed to obtain the Q field related to the path, as shown in FIG. Figure 5 The figure shows the path-related depth domain Q field contour map. The horizontal axis represents the CDP gather, the vertical axis represents the depth, and the contour lines inside the figure represent the offset. Figure 5 a)-f) in the figure are the Q fields corresponding to different offset groups, with offsets of 50m-783m, 783m-1516m, 1516m-2249m, 2249m-2982m, 2982m-3715m and 3715m-4450m, respectively. The path-dependent Q field distribution diagram is shown, and it can be seen that the Q fields at different offsets are non-uniformly distributed. Using the path-dependent Q field for viscoelastic prestack depth migration, such as Figure 6 The figure shows the viscoelastic prestack depth migration result based on the path-dependent equivalent Q field. Figure 4 The conventional prestack depth migration shown has a higher resolution. Further, the conventional prestack depth migration method and the viscoelastic prestack depth migration method of the present invention are used to perform viscoelastic prestack depth migration at different depths. Figure 7 The figure shows the results of viscoelastic prestack depth migration (1500m-4000m) of the 415th line of the imaging space in the work area. Figure 7 a) and b) are imaging results of the conventional method and the imaging results of the method of the present invention, respectively; Figure 8 The figure shows the results of viscoelastic prestack depth migration (2500m-5000m) of the 415th line of the imaging space in the work area. Figure 8 a) and b) are imaging results of the conventional method and the imaging results of the method of the present invention, respectively; Fig. 9 The following is the result of viscoelastic prestack depth migration (4200-6000 m) of the 415th line of imaging space in the work area. Fig. 9 a) and b) in the figure are the imaging results of the conventional method and the imaging results of the method of the present invention, respectively. Figure 7-Figure 9 The imaging results from shallow to deep layers shown in the figure show that the resolution of the imaging results obtained by the method of the present invention is significantly improved compared with the conventional method, and the phase axes are more continuous and well separated. Fig.10 a)-c) in the figure show respectively Figure 7-Figure 9 The spectrum comparison of the prestack depth migration results of the corresponding window shows that, compared with the conventional method, the effective bandwidth of the viscoelastic migration results of the method of the present invention has been significantly improved. The above simulation experiments fully demonstrate the effectiveness and feasibility of the method of the present invention.

[0118] like Fig.11As shown, corresponding to the above-mentioned viscoelastic prestack depth migration method based on path-dependent equivalent Q value, an embodiment of the present invention further provides a viscoelastic prestack depth migration system based on path-dependent equivalent Q value, and the above-mentioned viscoelastic prestack depth migration system based on path-dependent equivalent Q value includes:

[0119] An information acquisition module 111 is used to collect pre-stack seismic data of a target work area;

[0120] A constant Q scanning gather generation module 112 is used to perform inverse Q filtering on the pre-stack seismic data to obtain a depth domain constant Q scanning gather that varies with the offset group;

[0121] A constant Q value generating module 113 is used to perform a constant Q scan on the constant Q scan gather according to a preset time window, and determine a constant Q value of an offset group corresponding to each time window;

[0122] A path-dependent equivalent Q field construction module 114 is used to establish a path-dependent equivalent Q field based on the constant Q value of all the offset groups;

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

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

[0125] Based on the above embodiments, the present invention further provides an intelligent terminal, whose principle block diagram can be shown as follows: Fig.12 As shown. The above-mentioned intelligent terminal includes a processor, a memory, a network interface and a display screen connected through a system bus. Among them, the processor of the intelligent terminal is used to provide computing and control capabilities. The memory of the intelligent terminal includes a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system and a viscoelastic prestack depth migration program based on a path-related equivalent Q value. The internal memory provides an environment for the operation of the operating system and the viscoelastic prestack depth migration program based on a path-related equivalent Q value in the non-volatile storage medium. The network interface of the intelligent terminal is used to communicate with an external terminal through a network connection. When the viscoelastic prestack depth migration program based on a path-related equivalent Q value is executed by the processor, the steps of any one of the above-mentioned viscoelastic prestack depth migration methods based on a path-related equivalent Q value are implemented. 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 will understand that Fig.12The principle block diagram shown in the figure is only a block diagram of a partial structure related to the scheme of the present invention, and does not constitute a limitation on the smart terminal to which the scheme of the present invention is applied. The specific smart terminal may include more or fewer components than those shown in the figure, or combine certain components, or have a different arrangement of components.

[0127] In one embodiment, a smart terminal is provided. The smart terminal includes a memory, a processor, and a viscoelastic prestack depth migration program based on a path-correlated equivalent Q value stored in the memory and executable on the processor. When the viscoelastic prestack depth migration program based on a path-correlated equivalent Q value is executed by the processor, the steps of any one of the viscoelastic prestack depth migration methods based on a path-correlated equivalent Q value provided in an embodiment of the present invention are implemented.

[0128] An embodiment of the present invention further provides a computer-readable storage medium, on which is stored a viscoelastic prestack depth migration program based on a path-correlated equivalent Q value. When the viscoelastic prestack depth migration program based on a path-correlated equivalent Q value is executed by a processor, the steps of any one of the viscoelastic prestack depth migration methods based on a path-correlated equivalent Q value provided in an embodiment of the present invention are implemented.

[0129] It should be understood that the serial numbers of the steps in the above embodiments do not imply a sequence of execution, and the execution sequence of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of the present invention.

[0130] Those skilled in the art can clearly understand that, for the convenience and simplicity of description, only the division of the above-mentioned functional units and modules is used as an example for illustration. In practical applications, the above-mentioned function allocation can be completed by different functional units and modules as needed, that is, the internal structure of the above-mentioned device can be divided into different functional units or modules to complete all or part of the functions described above. The functional units and modules in the embodiment can be integrated in a processing unit, or each unit can exist physically separately, or two or more units can be integrated in one unit. The above-mentioned integrated unit can be implemented in the form of hardware or in the form of software functional units. In addition, the specific names of the functional units and modules are only for the convenience of distinguishing each other, and are not used to limit the scope of protection of the present invention. The specific working process of the units and modules in the above-mentioned system can refer to the corresponding process in the aforementioned method embodiment, which will not be repeated here.

[0131] In the above embodiments, the description of each embodiment has its own emphasis. For parts that are not described or recorded in detail in a certain embodiment, reference can be made to the relevant descriptions of other embodiments.

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

[0133] In the embodiments provided by the present invention, it should be understood that the disclosed apparatus / terminal equipment and method can be implemented in other ways. For example, the apparatus / terminal equipment embodiments described above are only illustrative, for example, the division of the above modules or units is only a logical function division, and in actual implementation, other division methods can be used, for example, multiple units or components can be combined or integrated into another system, or some features can be ignored or not executed.

[0134] The embodiments described above are only used to illustrate the technical solutions of the present invention, rather than to limit the same. Although the present invention has been described in detail with reference to the aforementioned embodiments, it should be understood by those skilled in the art that the technical solutions described in the aforementioned embodiments may still be modified, or some of the technical features therein may be replaced by equivalents. However, these modifications or replacements do not deviate from the spirit and scope of the technical solutions of the embodiments of the present invention, and should be included in the protection scope of the present invention.

Claims

1. A viscoelastic prestack depth migration method based on path-dependent equivalent Q value, characterized in that: The following steps are involved: Collect pre-stack seismic data in the target work area; Performing inverse Q filtering on the pre-stack seismic data to obtain a depth domain constant Q scanning gather that varies with the offset group; Performing a constant Q scan on the constant Q scan gathers according to a preset time window to determine a constant Q value of an offset group corresponding to each time window; establishing a path-dependent equivalent Q field based on a constant Q value for all offset groups; Viscoelastic prestack depth migration is performed based on the path-dependent equivalent Q field to obtain imaging results.

2. The viscoelastic prestack depth migration method based on path-dependent equivalent Q value according to claim 1, characterized in that: The step of performing a constant Q scan on the constant Q scan gathers according to a preset time window to determine a constant Q value of an offset group corresponding to each time window includes: Dividing the constant-Q scan gathers according to preset time windows, and determining a number of Q values ​​corresponding to each of the time windows; The superimposed frequency spectra corresponding to all the Q values ​​in each window are calculated, and based on the superimposed frequency spectra, the constant Q value of the offset group corresponding to each time window is determined.

3. The viscoelastic prestack depth migration method based on path-dependent equivalent Q value according to claim 1, characterized in that: The establishing of a path-dependent equivalent Q field based on the constant Q value of all the offset groups comprises: Based on the parabolic relationship between Q value and offset, a depth domain Q value model is established; Determining an optimal parameter combination of the depth domain Q value model corresponding to each time window based on the constant Q value of the offset group; Based on the optimal parameter combination, a path-dependent equivalent Q field is established.

4. The viscoelastic prestack depth migration method based on path-dependent equivalent Q value according to claim 3 is characterized in that: The determining, based on the constant Q value of the offset group, an optimal parameter combination of the depth domain Q value model corresponding to each time window comprises: Construct an objective function based on the constant Q spectrum and path-dependent Q spectrum corresponding to each time window; Based on the constant Q value of the offset group, a preset optimization algorithm is used to solve the optimal solution of the objective function corresponding to each time window; Based on the optimal solution, the best parameter combination of the depth-domain Q-value model corresponding to each time window is determined.

5. The viscoelastic prestack depth migration method based on path-dependent equivalent Q value according to claim 4 is characterized in that: The objective function is constructed based on the constant Q spectrum and the path-related Q spectrum corresponding to each time window, including: Calculate the center frequency, main low-frequency component, bandwidth and spectrum area of ​​the constant Q spectrum and path-related Q spectrum corresponding to each time window; Based on the center frequency, the main low-frequency component, the bandwidth and the spectrum area, constructing superimposed spectrum indicators and offset spectrum indicators corresponding to the constant Q spectrum and the path-related Q spectrum respectively; Based on all the superimposed spectrum indicators and the offset spectrum indicators, an objective function is constructed.

6. The viscoelastic prestack depth migration method based on path-dependent equivalent Q value according to claim 5, characterized in that: The method of solving the optimal solution of the objective function corresponding to each time window by using a preset optimization algorithm based on the constant Q value of the offset group includes: Initializing parameters in the depth-domain Q-value model to obtain initial values; Obtaining an initial solution of the objective function based on a superimposed spectrum and an offset spectrum corresponding to the constant-Q spectrum; Based on the initial value and the initial solution, a probability-based global optimization algorithm is used to find the optimal solution of the objective function corresponding to each time window.

7. The viscoelastic prestack depth migration method based on path-dependent equivalent Q value according to claim 3, characterized in that: The establishing of a path-dependent equivalent Q field based on the optimal parameter combination includes: Solving the depth domain Q value model based on the optimal parameter combination corresponding to each of the time windows to obtain the optimal Q value corresponding to each of the time windows; All the optimal Q values ​​are interpolated and smoothed to establish a path-dependent equivalent Q field.

8. A viscoelastic prestack depth migration system based on path-dependent equivalent Q value, characterized in that: The system comprises: Information acquisition module, used to collect pre-stack seismic data of the target work area; A constant Q scanning gather generation module is used to perform inverse Q filtering on the pre-stack seismic data to obtain a depth domain constant Q scanning gather that varies with the offset group; A constant Q value generating module, used for 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; A path-dependent equivalent Q field building module is used to build a path-dependent equivalent Q field based on the constant Q value of all the offset groups; The viscoelastic prestack depth migration module is used to perform viscoelastic prestack depth migration based on the path-related equivalent Q field to obtain imaging results.

9. An intelligent terminal, characterized in that: The intelligent terminal includes a memory, a processor, and a viscoelastic prestack depth migration program based on a path-correlated equivalent Q value stored in the memory and executable on the processor. When the viscoelastic prestack depth migration program based on a path-correlated equivalent Q value is executed by the processor, the steps of the viscoelastic prestack depth migration method based on a path-correlated equivalent Q value are implemented as described in any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores a viscoelastic prestack depth migration program based on a path-dependent equivalent Q value. When the viscoelastic prestack depth migration program based on a path-dependent equivalent Q value is executed by a processor, the steps of the viscoelastic prestack depth migration method based on a path-dependent equivalent Q value as described in any one of claims 1 to 7 are implemented.

Images