Velocity Spectrum Optimization Method for Seismic Data

Abstract

This application relates to a method for optimizing the velocity spectrum of seismic data, which includes obtaining the instantaneous amplitude of seismic data; performing empirical mode decomposition on the instantaneous amplitude, screening the decomposed intrinsic modes, selecting the intrinsic modes that are beneficial for expressing the velocity spectrum and superimposing them to generate new seismic data, and calculating the velocity spectrum using the new seismic data. This method obtains the instantaneous amplitude of seismic data based on the Hilbert transform, performs EMD decomposition on the instantaneous amplitude, screens the decomposed intrinsic modes, selects the intrinsic modes that are beneficial for expressing the velocity spectrum and superimposes them to generate new seismic data. The spectrum of this seismic data shows that the effective frequency band moves towards the low frequency. Calculating the velocity spectrum using the frequency-shifted seismic data can effectively improve the accuracy of the velocity spectrum. A high-precision velocity spectrum is beneficial to the accuracy of velocity analysis and reduces its multi-solution nature, thereby achieving the purpose of improving the imaging quality of seismic data.

Description

Technical Field

[0001] The present invention relates to the field of seismic data processing, and particularly to a method for optimizing the velocity spectrum of seismic data. Background Art

[0002] In the process of seismic data imaging, the accuracy and precision of formation velocity are key elements for imaging, and velocity analysis is the most commonly used means to obtain the velocity of underground reflection waves. Velocity analysis is based on quality control means such as velocity spectra, super gathers, small stacks, and stacked profiles as reference bases for velocity picking, among which the velocity spectrum is the most core reference basis in the velocity analysis process. Performing a velocity scan on the common midpoint prestack gathers, the change in the dynamic correction stack energy relative to the scanning velocity is called the velocity spectrum. The closer the scanning velocity corresponding to the focused position of the velocity spectrum energy cluster is to the true velocity at that position, the better the velocity spectrum has concentrated energy clusters longitudinally and obvious velocity trends, and is single and stable laterally. The main factors affecting the accuracy and resolution of the velocity spectrum include: (1) offset distribution; (2) fold number; (3) signal-to-noise ratio; (4) muting; (5) the bandwidth of seismic data, etc.

[0003] Marine seismic data are affected by rugged formations, marine multiples, and diffracted waves, etc., and the velocity spectrum energy clusters often show phenomena such as divergence and multiple energy clusters at the same reflection time, increasing the difficulty of velocity analysis. Improving the accuracy of the velocity spectrum can directly improve the accuracy of velocity modeling and reduce its non-uniqueness, and further improve the imaging quality of seismic data. Therefore, how to improve the accuracy of the velocity spectrum is a technical problem to be continuously solved. Summary of the Invention

[0004] The purpose of the present invention is to provide a method for optimizing the velocity spectrum of seismic data in view of the above deficiencies.

[0005] One technical solution adopted by the present invention to solve the technical problem is to provide a method for optimizing the velocity spectrum of seismic data, obtaining the instantaneous amplitude of seismic data; performing empirical mode decomposition on the instantaneous amplitude, screening the decomposed intrinsic modes, selecting the intrinsic modes that are beneficial to expressing the velocity spectrum and superimposing them to generate new seismic data, and calculating the velocity spectrum using the new seismic data;

[0006] The specific steps are as follows:

[0007] S1: Processing the seismic data using Hilbert transform and calculating the instantaneous amplitude spectrum;

[0008] S2: Performing empirical mode decomposition on the amplitude spectrum of the instantaneous amplitude;

[0009] S3: Screening the decomposed intrinsic modes, selecting the intrinsic modes that are beneficial to expressing the velocity spectrum and superimposing them to generate new seismic data;

[0010] S4: The effective frequency band shifts to lower frequencies in the new seismic data, and the velocity spectrum is calculated using the frequency-shifted seismic data.

[0011] Furthermore, in step S1, for a real-valued seismic signal x(t), the analytic signal c(t) is constructed using formula (1), the Hilbert transform signal H(t) is constructed using formula (2), the instantaneous amplitude A(t) is calculated using formula (3), and the instantaneous phase φ(t) is calculated using formula (4);

[0012] The calculation formulas are as follows:

[0013] Formula (1) is c(t) = x(t) + iH(t) = A(t)e iφ(t) , where i is the imaginary unit

[0014] Formula (2) is

[0015] Formula (3) is

[0016] Formula (4) is φ(t) = Arg[c(t)].

[0017] Furthermore, in step S2, first, the empirical mode decomposition is performed on the instantaneous amplitude signal of formula (4) using formula (5). The empirical mode decomposes the non-linear and non-stationary signal into several intrinsic mode functions and a residue;

[0018] Formula (5) is where IMF k (t) represents the intrinsic mode function; K represents the number of intrinsic mode functions; Res(t) represents the residue.

[0019] Furthermore, for the intrinsic mode function IMF k (t) within the entire time range, the number of local extreme points and zero-crossing points differs by at most one, and at any moment, the envelopes of the local maximum and local minimum must average to zero.

[0020] Furthermore, in step S3, according to the frequency band range of the actual seismic data and the frequency range of the intrinsic mode function IMF k (t), starting from the first intrinsic mode function, L consecutive intrinsic mode functions that are beneficial to the velocity spectrum are selected and superimposed using formula (6) to construct the new seismic data x new (t);

[0021] Formula (6) is where L ≤ K;

[0022] Subsequently, the new residue Res is obtained using formula (7)new (t);

[0023] Equation (7) is Res new (t) = A(t) - x new (t).

[0024] Furthermore, in step S4, within the common midpoint domain, the newly constructed seismic data x new (t) is used to calculate the stacking velocity spectrum. The specific calculation method is as follows:

[0025] S4.1: Take out several adjacent CMP gathers near the velocity point to be analyzed;

[0026] S4.2: Average the taken adjacent CMP gathers to form a super gather;

[0027] S4.3: Select several constant velocities to perform NMO correction and stacking on the super gather. One constant velocity corresponds to one stacked section. The stacked sections obtained from different velocities constitute the CVS section;

[0028] S4.4: Interpret the ideal stacking velocity from the CVS image using the maximum stacking energy criterion. The velocity corresponding to the strongest post-stack event amplitude and good continuity is the required stacking velocity;

[0029] S4.5: Modify the value of L, and repeat steps S3 and S4 until at 80% of the vertical time of the velocity spectrum and there is only a single energy cluster horizontally, and the width of the energy cluster is less than 300 m / s, then output the optimized velocity spectrum.

[0030] Furthermore, in step S4.4, when the velocity is appropriate, after NMO correction, the reflection event is flattened. At this time, the stacking effect is the best, the stacking amplitude shows a maximum value, and the stacking energy is the strongest;

[0031] The stacking energy calculation formula (8) is as follows: In the formula, M is the number of horizontal seismic traces in the CMP gather; N is the number of samples within the vertical time window for calculating the stacking velocity; x' new (t m,n ) is the amplitude value of the nth sample of the mth trace after NMO correction of x new (t); is the average amplitude within the time window, that is, the stacking energy.

[0032] Furthermore, before step S1, the seismic gather data needs to be input first.

[0033] Another technical solution adopted by the present invention to solve the technical problem is to provide a terminal, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, the above-mentioned velocity spectrum optimization method for seismic data is realized.

[0034] Compared with the prior art, the present invention has the following beneficial effects: The instantaneous amplitude of seismic data is obtained based on the Hilbert transform, the empirical mode decomposition (EMD) is performed on the instantaneous amplitude, the intrinsic mode functions (IMFs) after decomposition are screened, and the IMFs beneficial to expressing the velocity spectrum are selected and superimposed to generate new seismic data. The spectrum of this new seismic data shows that the effective frequency band shifts to the low frequency. Calculating the velocity spectrum using the frequency-shifted seismic data can effectively improve the accuracy of the velocity spectrum; the high-precision velocity spectrum is beneficial to the accuracy of velocity analysis and reduces its multi-solution property, thereby achieving the purpose of improving the imaging quality of seismic data. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] The following further describes the present invention in conjunction with the drawings.

[0036] Figure 1 It is a flow chart of the method.

[0037] Figure 2 It is the original stacked data and the instantaneous amplitude.

[0038] Figure 3 It is the EMD decomposition diagram of the instantaneous amplitude of the stacked data. (a) is the instantaneous amplitude; (b) is the intrinsic mode function IMF 1 (t); (c) is the intrinsic mode function IMF 2 (t); (d) is the intrinsic mode function IMF 3 (t); (e) is the intrinsic mode function IMF 4 (t); (f) is the intrinsic mode function IMF 5 (t); (g) is the residual Res(t) (amplitude amplified by 10 7 ).

[0039] Figure 4 It is the stacked data and the residual after frequency shift processing. On the left, the intrinsic mode functions IMF 1 (t) and IMF 2 (t) are added to construct the new seismic data x new (t); on the right is the new residual Res new (t).

[0040] Figure 5 It is the comparison diagram of the spectra of the stacked data before and after frequency shift processing.

[0041] Figure 6 It is the comparison of the common midpoint gathers before and after frequency shift processing.

[0042] Figure 7 It is a comparison diagram of the spectra of the common midpoint gather before and after frequency shift processing.

[0043] Figure 8 It is a comparison diagram of the velocity spectra before and after frequency shift processing. Specific implementation manners

[0044] The preferred implementation manners of the present invention will be described in more detail below. However, it should be understood that the present invention can be implemented in various forms and should not be limited by the implementation manners described herein. On the contrary, these implementation manners are provided to make the present invention more thorough and complete, and to fully convey the scope of the present invention to those skilled in the art. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0045] Embodiment 1:

[0046] The Hilbert transform is a commonly used method in seismic signal processing. Seismic signals are real signals, and their corresponding complex signals need to be calculated in the complex signal analysis method. Taking the seismic signal to be analyzed as the real part and the Hilbert transform of the seismic signal to be analyzed as the imaginary part to form a complex signal, which can be used to calculate the instantaneous attributes of the signal. The three instantaneous attributes (instantaneous amplitude, instantaneous frequency, instantaneous phase) of seismic data obtained based on the Hilbert transform can be used for the frequency study of seismic data.

[0047] The velocity spectrum is low-frequency relative to seismic data. Then, improving the low-frequency components of the seismic data used to calculate the velocity spectrum can improve the quality of the velocity spectrum.

[0048] Based on the above signal analysis technical means, Embodiment 1 provides a method for optimizing the velocity spectrum of seismic data based on the Hilbert transform to improve the effective bandwidth of the common midpoint gather, especially to strengthen the low-frequency components of the gather, so as to achieve the purpose of optimizing the velocity spectrum.

[0049] As Figure 1 shown, a method for optimizing the velocity spectrum of seismic data includes obtaining the instantaneous amplitude of seismic data; performing empirical mode decomposition on the instantaneous amplitude, screening the decomposed eigenmodes, selecting the eigenmodes that are beneficial to expressing the velocity spectrum and superimposing them to generate new seismic data, and calculating the velocity spectrum using the new seismic data;

[0050] The specific steps are as follows:

[0051] S1: Process the seismic data using the Hilbert transform and calculate the instantaneous amplitude spectrum;

[0052] In step S1, for a real-valued seismic signal x(t), an analytic signal c(t) is constructed using formula (1), a Hilbert transform signal H(t) is constructed using formula (2), an instantaneous amplitude A(t) is calculated using formula (3), and an instantaneous phase φ(t) is calculated using formula (4);

[0053] The calculation formulas are as follows:

[0054] Formula (1) is c(t) = x(t) + iH(t) = A(t)e iφ(t) , where i is the imaginary unit

[0055] Formula (2) is

[0056] Formula (3) is

[0057] Formula (4) is φ(t) = Arg[c(t)];

[0058] S2: Perform empirical mode decomposition on the amplitude spectrum of the instantaneous amplitude;

[0059] In step S2, first, the instantaneous amplitude signal of formula (4) is subjected to empirical mode decomposition using formula (5). Empirical mode decomposes a non-linear and non-stationary signal into several intrinsic mode functions and a residue;

[0060] Formula (5) is where IMF k (t) represents the intrinsic mode function; K represents the number of intrinsic mode functions; Res(t) represents the residue;

[0061] The intrinsic mode function IMF k (t) has at most a difference of one between the number of local extreme points and zero-crossing points within the entire time range, and at any time point, the average of the envelopes of the local maximum (upper envelope) and the local minimum (lower envelope) must be zero;

[0062] S3: Screen the decomposed intrinsic mode functions, select the intrinsic mode functions that are beneficial to expressing the velocity spectrum, and superimpose them to generate new seismic data;

[0063] In step S3, according to the frequency band range of the actual seismic data and the frequency range of the intrinsic mode function IMF k (t), starting from the first intrinsic mode function, L consecutive intrinsic mode functions that are beneficial to the velocity spectrum are selected (this step aims to optimize the velocity spectrum and must start from the first intrinsic mode function without interruption. For example, retain the first 5 or the first 3 mode functions), and formula (6) is used for superposition to construct new seismic data x new (t);

[0064] Formula (6) is where L ≤ K;

[0065] Subsequently, a new residual Res new (t) is obtained using formula (7);

[0066] Formula (7) is Res new (t) = A(t) - x new (t);

[0067] Compared with the amplitude spectrum of the original seismic data x(t), the effective frequency of the new seismic data x new (t) shifts towards the low-frequency end, and the scale of the shift depends on the value of L. The whole process of the above steps is called frequency shift processing;

[0068] S4: It is shown that the effective frequency band shifts towards the low frequency on the new seismic data, and the velocity spectrum is calculated using the seismic data after frequency shift;

[0069] In step S4, within the common midpoint domain, the newly constructed seismic data x new (t) is used to calculate the stacking velocity spectrum, and the specific calculation method is as follows:

[0070] S4.1: Take out several adjacent CMP gathers near the velocity point to be analyzed;

[0071] In practical applications, usually 5 or 7 adjacent CMP gathers are selected;

[0072] S4.2: The adjacent CMP gathers taken out are averaged to form a super gather;

[0073] S4.3: Select several constant velocities (for example: 1500 m / s - 6000 m / s, at intervals of 100 m / s) to perform dynamic correction and stacking on the super gather. One constant velocity corresponds to one stacked section, and the stacked sections obtained from different velocities constitute the CVS section; the stacking energy of the CVS section increases or decreases with the change of the several constant velocities, and the change of the stacking energy with respect to the velocity is the velocity spectrum;

[0074] S4.4: Interpret the ideal stacking velocity from the CVS image using the maximum stacking energy criterion. The velocity corresponding to the strongest amplitude of the post-stack event axis and good continuity is the required stacking velocity; the amplitude of the stacked record changes with the stacking velocity, and this is the stacking velocity spectrum;

[0075] In step S4.4, when the velocity is appropriate, after dynamic correction, the reflection wave event axis is flattened, and at this time, the stacking effect is the best, the stacking amplitude appears as a maximum value, and the stacking energy is the strongest;

[0076] The stacking energy calculation formula (8) is as follows: In the formula, M is the number of lateral seismic traces in the CMP gather; N is the number of samples in the longitudinal time window for calculating the stacking velocity; x' new (t m,n ) is the amplitude value of the nth sample of the mth trace after NMO correction of x new (t); is the average amplitude within the time window, i.e., the stacking energy;

[0077] S4.5: Modify the value of L, and repeat steps S3 and S4 until at 80% of the longitudinal time of the velocity spectrum, there is only a single energy cluster in the transverse direction, and the width of the energy cluster is less than 300 m / s, then output the optimized velocity spectrum.

[0078] In this embodiment, before step S1, it is necessary to input the seismic gather data first.

[0079] This method obtains the instantaneous amplitude of seismic data based on the Hilbert transform, performs EMD decomposition on the instantaneous amplitude, screens the decomposed intrinsic mode functions, selects the intrinsic mode functions beneficial to expressing the velocity spectrum and stacks them to form new seismic data. The spectrum of this seismic data shows that the effective frequency band shifts to the low frequency. Calculating the velocity spectrum using the frequency-shifted seismic data can effectively improve the accuracy of the velocity spectrum; a high-precision velocity spectrum is beneficial to the accuracy of velocity analysis and reduces its ambiguity, thereby achieving the purpose of improving the imaging quality of seismic data.

[0080] The following provides two verification examples of Example 1:

[0081] Verification Example 1:

[0082] Taking the stacked data of an actual marine streamer seismic data as an example, the recording length of this data is 7 s, the sampling rate is 1 ms, the trace interval is 3.125 m, and the effective frequency band range is about 4 - 120 Hz. Generally, the longitudinal and transverse resolutions are both good. This method is used to perform frequency shift processing on this data. Figure 2 is the original stacked data and the instantaneous amplitude diagram. Figure 3 is the EMD decomposition diagram of the instantaneous amplitude. It can be seen that the first decomposed intrinsic mode function IMF 1 (t) has the highest frequency, and the frequencies of the subsequent EMD-decomposed intrinsic mode functions gradually decrease. The frequency of the intrinsic mode function IMF 5 (t) decomposed in the 5th time is less than 4 Hz, and the amplitude of the residual Res(t) is very small and can only be seen by the naked eye after being amplified by 10 7 times. Add the intrinsic mode functions IMF 1 (t) and IMF 2 (t) to construct new seismic data x new (t) and calculate the new residual Res new (t). Figure 4are the stacked data and residuals before and after frequency shift processing. Figure 5 is the spectrum comparison chart before and after frequency shift processing. The stacked data after frequency shift processing shows obvious low-frequency energy. Combining with Figure 5 the shown spectrum comparison chart, the spectrum of the instantaneous amplitude presents an exponential form, with strong low-frequency in the spectrum; after frequency shift processing, the effective frequency band moves towards the low-frequency end, which is consistent with the obvious low-frequency energy characteristics shown in the stacked data after processing.

[0083] As can be seen from the above, Verification Example 1 proves that frequency shift processing can effectively move the frequency band of seismic data towards the low-frequency end.

[0084] Verification Example 2:

[0085] Taking the prestack common midpoint data of an ocean towed streamer in a rugged formation as an example, the fold is 80 times, the record length is 12s, the sampling rate is 2ms, the trace interval is 12.5m, and the effective frequency band is relatively narrow (about 4 - 90Hz). Figure 6 is the comparison chart of common midpoint gathers before and after frequency shift processing. Before processing, there are many effective reflection layers but lack wave group characteristics, showing the characteristics of narrow-band signals; after processing, the contrast between strong and weak amplitudes is obvious, and the wave group relationship is clear and definite. Figure 7 is the spectrum comparison chart of common midpoint gathers before and after frequency shift processing. The effective frequency band moves towards the low-frequency end, effectively increasing the low-frequency information of the gathers. Figure 8 is the Figure 6 stacking velocity spectrum calculated from the gather data of. The effect is significant before and after the optimization of the velocity spectrum. Before optimization, affected by the rugged formation, the energy clusters are divergent, there are multiple energy clusters at the same reflection time, and there is multi-solution in velocity picking, increasing the velocity error. After optimization, the energy clusters are focused, the velocity analysis is simple and the accuracy is greatly improved, laying a solid foundation for fine velocity modeling.

[0086] Example 2:

[0087] Example 2 provides a terminal, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the above-mentioned method for optimizing the velocity spectrum of seismic data.

[0088] Unless otherwise specifically stated, the relative arrangements of components and steps, numerical expressions, and numerical values set forth in these embodiments do not limit the scope of the present application. In all the examples shown and discussed here, any specific value should be construed as merely exemplary and not as a limitation. Therefore, other examples of the exemplary embodiments may have different values. It should be noted that: like reference numerals and letters denote like items, so once an item is defined, it does not need to be further discussed subsequently.

[0089] In the description of the present application, it should be understood that the orientation or positional relationship indicated by orientation words such as "front, rear, upper, lower, left, right", "lateral, vertical, perpendicular, horizontal" and "top, bottom", etc. is usually based on the shown orientation or positional relationship, and is only for the convenience of describing the present application and simplifying the description. Without contrary statements, these orientation words do not indicate or imply that the device or element referred to must have a specific orientation or be constructed and operated in a specific orientation, so it cannot be understood as a limitation on the protection scope of the present application; the orientation words "inside, outside" refer to the inside and outside relative to the contour of each component itself.

[0090] For the convenience of description, spatial relative terms such as "above...", "over...", "on the upper surface of...", "upper...", etc. can be used here to describe the spatial positional relationship of features. It should be understood that spatial relative terms are intended to include different orientations in use or operation in addition to the described orientation.

[0091] In addition, it should be noted that the use of words such as "first", "second", etc. to limit components is only for the convenience of distinguishing the corresponding components. Without additional statements, the above words have no special meaning, so it cannot be understood as a limitation on the protection scope of the present application.

[0092] If the present application discloses or involves components or structural members that are fixedly connected to each other, then, unless otherwise stated, the fixed connection can be understood as: a detachable fixed connection (for example, connected by bolts or screws), or can also be understood as: a non-detachable fixed connection (for example, riveting, welding). Of course, the mutually fixed connection can also be replaced by an integral structure (for example, manufactured by integral forming using a casting process) (except when it is obviously impossible to adopt the integral forming process).

[0093] The above preferred embodiments have further elaborated on the purpose, technical solution and advantages of the present invention. It should be understood that the above are only the preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

Claims

1. A method for optimizing the velocity spectrum of seismic data, characterized in that: obtain the instantaneous amplitude of the seismic data; perform empirical mode decomposition on the instantaneous amplitude, screen the decomposed intrinsic mode functions, select the intrinsic mode functions that are beneficial to expressing the velocity spectrum and superimpose them to generate new seismic data, and calculate the velocity spectrum using the new seismic data; The specific steps are as follows: S1: Process the seismic data using the Hilbert transform and calculate the instantaneous amplitude spectrum; S2: Perform empirical mode decomposition on the amplitude spectrum of the instantaneous amplitude; S3: Screen the decomposed intrinsic mode functions, select the intrinsic mode functions that are beneficial to expressing the velocity spectrum and superimpose them to generate new seismic data; S4: It is shown that the effective frequency band shifts to the low frequency on the new seismic data, and calculate the velocity spectrum using the frequency-shifted seismic data.

2. The method for optimizing the velocity spectrum of seismic data according to claim 1, characterized in that: In step S1, for a real-valued seismic signal , an analytic signal is constructed using formula (1) , a Hilbert transform signal is constructed using formula (2) , the instantaneous amplitude is calculated using formula (3) , the instantaneous phase is calculated using formula (4) : The calculation formula is as follows: Formula (1) is , where is the imaginary unit ; Formula (2) is ; Formula (3) is ; Formula (4) is .

3. The method for optimizing the velocity spectrum of seismic data according to claim 2, characterized in that: In step S2, first, perform empirical mode decomposition on the instantaneous amplitude signal of formula (4) using formula (5), and the empirical mode decomposes the non-linear and non-stationary signal into several intrinsic mode components and a residue; Equation (5) is , where represents the intrinsic mode component; K represents the number of intrinsic mode components; represents the residual.

4. The method for optimizing the velocity spectrum of seismic data according to claim 3, characterized in that: Intrinsic mode component Over the entire time range, the number of local extreme points and zero-crossing points differs by at most one, and at any given time point, the envelopes of the local maxima and local minima must average to zero.

5. The method for optimizing the velocity spectrum of seismic data according to claim 3, characterized in that: In step S3, according to the frequency band range of the actual seismic data and the frequency range of the intrinsic mode components , starting from the first intrinsic mode component, L consecutive intrinsic mode components beneficial to the velocity spectrum are screened out, and the new seismic data is constructed by superposition using formula (6) ; subsequently, the new residual is obtained using formula (7) ; Formula (6) is , where ; Formula (7) is .

6. The method for optimizing the velocity spectrum of seismic data according to claim 4, characterized in that: In step S4, within the common center point domain, the newly constructed seismic data is used to calculate the stacking velocity spectrum. The specific calculation method is as follows: S4.1: Take out several adjacent CMP gathers near the velocity point to be analyzed; S4.2: Average the taken adjacent CMP gathers to form a super gather; S4.3: Select several constant velocities to perform NMO correction and stacking on the super gather. One constant velocity corresponds to one stacked section, and the stacked sections obtained from different velocities constitute a CVS section; S4.4: Interpret the ideal stacking velocity from the CVS image using the maximum stacking energy criterion. The velocity corresponding to the strongest amplitude of the post-stack event and good continuity is the required stacking velocity; S4.5: Modify the value of L, and repeat steps S3 and S4 until at 80% of the vertical time of the velocity spectrum and only a single energy cluster exists in the horizontal direction, and the width of the energy cluster is less than 300 m / s, then output the optimized velocity spectrum.

7. The method for optimizing the velocity spectrum of seismic data according to claim 6, characterized in that: In step S4.4, when the velocity is appropriate, after NMO correction, the reflection event axis is flattened, and at this time the stacking effect is the best, the stacking amplitude appears as a maximum value, and the stacking energy is the strongest; The superposition energy calculation formula (8) is as follows: In the formula, M is the number of lateral seismic traces in the CMP gather; N is the number of samples within the longitudinal time window for calculating the stacking velocity; is the amplitude value of the nth sample of the mth trace after NMO correction; is the average amplitude within the time window, i.e., the superposition energy.

8. The method for optimizing the velocity spectrum of seismic data according to any one of claims 1-7, characterized in that: Before step S1, it is necessary to input the seismic gather data first.

9. A terminal, including a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: when the processor executes the computer program, it implements the method for optimizing the velocity spectrum of seismic data according to any one of claims 1-7.

Images