A method for identifying marine hydrates without typical seismic reflection characteristics

Abstract

The present invention discloses a method for identifying marine hydrates without typical seismic reflection characteristics, which includes steps such as optimizing pre-stack gathers, picking horizons, analyzing coherent velocity along horizons, extracting velocity characteristics, and interpreting velocity anomalies. Combining with the characteristics of seismic data of marine hydrates without typical seismic reflection characteristics, a technical process suitable for identifying such hydrates is designed. Through optimizing pre-stack gathers, picking horizons, and inverting coherent velocity, relatively accurate velocity information of the strata possibly containing hydrates can be obtained, and then through processing processes such as extracting velocity characteristics and interpreting velocity anomalies, the identification of marine hydrates without typical seismic reflection characteristics is realized. The present invention makes full use of seismic velocity information, is easy to implement, and provides an efficient and economical method for identifying marine hydrates without typical seismic reflection characteristics.

Description

Technical Field

[0001] The implementation of the present invention belongs to the technical field of marine natural gas hydrate resource exploration, and particularly relates to a method for identifying marine natural gas hydrates without typical seismic reflection characteristics. Background Art

[0002] Natural gas hydrate is a potential high-efficiency clean energy source. Geophysical methods are one of the effective methods for marine natural gas hydrate exploration and evaluation. The presence of natural gas hydrate makes the density of the sedimentary formation lower and the acoustic velocity higher, forming typical seismic reflection characteristics such as bottom-simulating reflector (BSR), polarity reversal, and amplitude blanking zone on seismic profiles. BSR is commonly found in submarine sedimentary formations and marine slopes. However, in some areas where hydrates are discovered, such as flat seabeds, BSR characteristics may not necessarily be shown on seismic profiles. The appearance of typical seismic reflection characteristics is related to factors such as the accumulation mode, occurrence state, occurrence depth, and free gas saturation of hydrates. How to identify hydrates in areas without typical seismic reflection characteristics such as BSR is a difficult problem.

[0003] When high-saturation natural gas hydrates are present in the submarine formation, the longitudinal wave velocity of the formation suddenly increases. If there is free gas in the underlying formation of the hydrate formation, the wave velocity will decrease or even reverse. The hydrate formation usually appears as a high-velocity geological body in a low-velocity background, and this geological body is approximately parallel to the seabed horizontally. For hydrate formations without typical reflection characteristics such as BSR, or thin hydrate layers immediately below the seabed, the velocity change characteristics contained in seismic data can be used to identify the hydrate formation. However, the occurrence states of marine hydrates may be thin-layered, massive, dispersed, veined, etc., which are diverse and complex. Only under suitable conditions can obvious velocity anomalies be seen in conventional velocity analysis.

[0004] Seismic velocity analysis is greatly affected by humans, with limited accuracy, and artificial velocity analysis is usually carried out at regular CDPs. At the same time, for the case where the hydrate layer is thin-layered and horizontally discontinuous, there is no obvious reflection interface between the hydrate layer and the surrounding formations, and velocity anomalies are difficult to be directly reflected in seismic reflection data. However, the seismic reflection characteristics represent the superposition response of adjacent reflected waves. Affected by the high-velocity hydrate formation, the velocity at the actual position calculated according to the coherence algorithm will show a local increase, and the increase amplitude is related to the velocity of the hydrate layer itself and its depth from the actual position.

[0005] Based on the above considerations, there is an urgent need to propose a method for identifying hydrates using seismic velocity information, combined with horizon picking and fine velocity analysis along the horizon method, to solve the problem of identifying marine natural gas hydrates without typical seismic reflection characteristics. Summary of the Invention

[0006] In view of the problem of identifying marine natural gas hydrates without typical seismic reflection characteristics, the present invention proposes a method for identifying marine natural gas hydrates without typical seismic reflection characteristics, which solves the problem of identifying hydrates without typical reflection characteristics based on seismic velocity information.

[0007] The present invention is implemented by the following technical solutions: A method for identifying marine hydrates without typical seismic reflection characteristics, comprising the following steps:

[0008] Step A: Optimize the pre-stack CMP gather;

[0009] Step B: Horizon picking: On the stacked section, pick up the reflection times of the sea floor and the top and bottom interfaces of the target horizon according to the seismic reflection characteristics;

[0010] Step C: Along-layer coherence velocity analysis: On the optimized pre-stack CMP gather, with the sea floor and horizon information picked up in Step B as constraints, obtain the velocity of the target horizon based on the coherence velocity inversion algorithm;

[0011] Step D: Velocity feature extraction: Combine the velocity difference between the target horizon and the overlying strata calculated in Step C, represent all negative values as velocity anomalies and positive values as normal velocities, and complete the calculation of the velocity differences for all survey lines;

[0012] Step E: Velocity anomaly interpretation: Combine the plane distribution map of the velocity differences of the target horizon in the work area obtained in Step D, and the velocity anomaly area corresponds to the hydrate distribution area.

[0013] Further, the specific implementation of Step B is as follows:

[0014] (1) On the stacked section, first automatically pick up the sea floor reflection time according to the energy ratio method. For the case where the direct wave and the sea floor reflection energy overlap, or the sea floor isochron is deformed due to noise, combine the energy ratio method and the correlation method to determine the sea floor reflection time;

[0015] (2) Based on whether it is parallel to the sea floor reflection isochron, and whether there is partial amplitude enhancement or weakening of the isochron, or the phenomenon of the isochron convexity upwards, determine the top and bottom interfaces of the target horizon.

[0016] Further, the along-layer coherence velocity inversion in Step C is specifically implemented as follows:

[0017] Step C1: According to the Dix formula, calculate the layer velocity using the seismic stacking velocity:

[0018]

[0019] where V nis the interval velocity between the nth and (n - 1)th interfaces, v n is the stacking velocity of the nth interface, T n is the zero-offset travel time of the nth interface;

[0020] Step C2: For a specific horizon, set a series of trial velocity values based on the initial velocity and at a certain velocity increment;

[0021] Step C3: For a given trial velocity value, calculate the depth of the interface by combining the reflection times of the top and bottom interfaces of the target horizon, and calculate the travel time of the reflected wave in the CMP gather through ray tracing;

[0022] Step C4: Calculate the coherence coefficient of the CMP gather along the calculated travel time curve;

[0023]

[0024] In the formula, S l,m represents the coherence coefficient of the lth CMP gather, τ represents the time window width, N represents the number of seismic traces in the gather, U(t i,m +j,x i ) represents the amplitude, x i is the offset of the ith trace, time (V n ) m represents the current trial velocity;

[0025] Step C5: Calculate the coherence coefficient corresponding to each trial velocity; the trial velocity corresponding to the maximum coherence coefficient is the interval velocity at this CMP position of the target horizon;

[0026] Step C6: Process multiple CMP gathers for velocity picking according to steps C1 - C5, obtain the interval velocities at multiple CMP positions, interpolate them to get the interval velocities at all CMP points, then perform smoothing processing to establish the initial model of the lateral variation of the interval velocity;

[0027] Step C7: Calculate the coherence function of all CMP gathers according to the current initial velocity model:

[0028]

[0029] In the formula, S T is the coherence function, T represents the number of CMP gathers.

[0030] Step C8: Continuously iteratively adjust the velocity model. When the error between the coherence function calculated in this iteration and the coherence function calculated in the previous iteration is less than the error threshold, the velocity in this iteration is the result of the coherence velocity inversion, and thus the velocity of the target horizon is obtained.

[0031] Further, the speed feature extraction in step D is specifically performed in the following manner:

[0032] If there is a thin layer of high-speed hydrate immediately below the seabed, the seabed speed shows a high-speed anomaly at the corresponding position, and the lateral distribution range of the high-speed anomaly indicates the lateral distribution range of the hydrate high-speed layer; for the target layer below the seabed, calculate the speed difference between the target layer and the overlying formation, and represent all negative values as speed anomalies and positive values as normal speeds.

[0033] Further, the optimization operations in step A include noise suppression, flattening of the gather, and excision of far-offset stretching.

[0034] Compared with the prior art, the advantages and positive effects of the present invention are as follows:

[0035] Based on the characteristics of marine hydrate seismic data without typical seismic reflection features, this solution designs a technical process suitable for the identification of such hydrates. Through processing procedures such as pre-stack gather optimization, horizon picking, coherent velocity inversion, speed anomaly calculation, and hydrate identification, it solves the problem of identifying marine hydrates without typical seismic reflection features, makes full use of seismic velocity information, obtains a relatively accurate speed of the possible hydrate-bearing layer, and conducts the identification of marine hydrates without typical seismic reflection features based on this, which is economical, efficient, convenient, and fast. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] Figure 1 It is a schematic diagram of the seismic stacking profile of the embodiment of the present invention, where S1 and S2 are the picked formation interfaces;

[0037] Figure 2 It is a schematic diagram of the speed value of the overlying formation of the target layer (the bottom interface is labeled as S1 in Figure 1 );

[0038] Figure 3 It is a schematic diagram of the speed value of the target layer (the top and bottom interfaces are labeled as S1 and S2 respectively in Figure 1 );

[0039] Figure 4 It is a schematic diagram of the speed value curves of the target layer and its overlying formation;

[0040] Figure 5 It is a schematic diagram of the speed difference curves of the target layer and its overlying formation;

[0041] Figure 6 It is a plan view of the speed difference between the bottom and top interfaces of the possible hydrate-bearing layer in the work area of the embodiment of the present invention, where the dark area is the possible location of hydrates. DETAILED DESCRIPTION OF THE INVENTION

[0042] To better understand the above objects, features, and advantages of the present invention, the present invention will be further described below in conjunction with the accompanying drawings and embodiments. In the following description, many specific details are set forth to fully understand the present invention. However, the present invention can also be implemented in other ways different from those described herein. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0043] This embodiment proposes a method for identifying marine natural gas hydrates without typical seismic reflection characteristics. Based on the seismic data after high-precision migration imaging, the specific processing steps are as follows:

[0044] Step A: Optimize the pre-stack CMP gather.

[0045] Step B: Horizon picking. On the stacked section, pick the reflection times of the seabed and the top and bottom interfaces of the target horizon according to the seismic reflection characteristics.

[0046] Step C: Along-layer coherence velocity analysis. On the optimized pre-stack CMP gather, with the information of the picked seabed and target horizon as constraints, obtain the velocity of the target horizon based on the coherence velocity inversion algorithm.

[0047] Step D: Velocity feature extraction. Calculate the velocity difference between the target horizon and the overlying strata. Represent all negative values as velocity anomalies and positive values as normal velocities; and complete the calculation of the velocity differences for all survey lines.

[0048] Step E: Velocity anomaly interpretation. Obtain the plane distribution map of the velocity differences of the target horizon in the work area. The velocity anomaly area corresponds to the hydrate distribution area.

[0049] Specifically, the solution of the present invention will be described in detail below:

[0050] In step A, when optimizing the pre-stack gather, it includes but is not limited to operations such as removing post-migration noise, flattening the gather, and stretching and cutting of far offset distances.

[0051] In step B, when performing horizon picking, the following specific method is adopted:

[0052] (1) On the stacked section, first automatically pick the seabed reflection time according to the energy ratio method:

[0053]

[0054] where E(a) is the energy ratio of the time window after time a to the time window before time a; U(i) is the amplitude at the i-th time, and L represents the time window length.

[0055] (2) For the case where the direct wave and the seabed reflection energy overlap, or the seabed event is distorted due to noise, it is necessary to combine the energy ratio method and the correlation method to determine the seabed reflection time. First, use the energy ratio method to pick up the approximate seabed reflection time. Taking this as a reference, calculate the correlation coefficient of the data to further determine the accurate seabed time:

[0056]

[0057] In the formula, R(a) is the similarity coefficient, τ is the time delay, and t1 and t2 are the start time and end time of correlation calculation respectively.

[0058] The principle of inline coherence velocity analysis in step C is as follows:

[0059] (1) According to Dix formula, use seismic stacking velocity to calculate interval velocity:

[0060]

[0061] In the formula, V n is the interval velocity between the nth and (n - 1)th reflection interfaces, v n is the stacking velocity of the nth reflection interface, and T n is the zero-offset travel time of the nth reflection interface.

[0062] (2) Take the velocity of the target horizon on the lth CMP gather as the initial velocity (V n )0;

[0063] (3) For the target horizon, taking the initial velocity (V n )0 as the reference, set 2L + 1 test velocity values (V n ) m =(V n )0 + mΔV, where m = -L, -L + 1,..., -1, 0, 1,..., L - 1, L;

[0064] (4) For a given test velocity value (V n ) m , combined with the picked reflection time T0 of the formation interface, calculate the depth of the interface, and calculate the reflection wave travel time (T0) m of the CMP gather through ray tracing;

[0065] (5) Using the test velocity (V n ) m and the travel time (T0) m , calculate the forward gather and calculate its coherence coefficient with the actual CMP gather:

[0066]

[0067] In the formula, S l,m represents the coherence coefficient of the m-th test velocity, τ represents the time window width, N represents the number of channels, U(t i,m +j,x i ) represents the amplitude, and x i is the offset of the i-th channel, and the time

[0068] (6) Calculate the coherence coefficients corresponding to 2L + 1 test velocities respectively, and take the test velocity corresponding to the maximum coherence coefficient as the layer velocity at this CMP position of the target layer.

[0069] (7) Process the CMP gather for velocity picking according to steps (1)-(6) to obtain the layer velocities at these CMP positions, interpolate them to obtain the layer velocities at all CMP points, perform smoothing processing, and establish an initial model of the lateral variation of the layer velocity;

[0070] (8) According to the current initial velocity model, calculate the coherence functions of all CMP gathers:

[0071]

[0072] In the formula, S T is the coherence function, and T represents the number of CMP gathers.

[0073] (9) Continuously perform iterative adjustment on the velocity model. When the error between the coherence function calculated in this iteration (the k-th iteration) and the coherence function calculated in the previous iteration (the k - 1-th iteration) is less than the error threshold, the velocity of the k-th iteration is the result of the coherence velocity inversion.

[0074] |S T,k -S T,k-1 |<ε(6)

[0075] In the formula, S T,k , S T,k-1 are the coherence function values calculated after the k-th and k - 1-th iterations respectively, and ε is the set error threshold, which is specifically set according to the actual situation or experience.

[0076] When performing velocity feature extraction in step D, the following method is adopted:

[0077] (1) If there is a high-speed hydrate thin layer immediately below the seabed, the seabed velocity shows a high-speed anomaly at the corresponding position, and the lateral distribution range of the high-speed anomaly indicates the lateral distribution range of the hydrate high-speed layer;

[0078] (2) For the target layer beneath the seabed, calculate the velocity difference between the target layer and the overlying strata. Represent all negative values as velocity anomalies (inverted), and positive values as normal velocities.

[0079] When interpreting velocity anomalies in step E, according to the planar distribution map of the velocity differences of the potentially hydrate-bearing layers in the work area obtained in step D, if the distribution of the velocity anomaly areas is spatially continuous to a certain extent, it indicates a relatively high reliability of this result; some discontinuous or sporadically distributed anomalies may be caused by discontinuous reflection event axes and need to be comprehensively judged by further combining seismic profiles and the geological conditions of the work area.

[0080] Specifically, as Figures 1 - 6 shown, Figure 1 in, on the seismic stack profile obtained after optimization processing, pick up the top interface of the target layer ( Figure 1 marked as S1 in Figure 1 ) and the bottom interface ( Figure 2 marked as S2 in Figure 3 ) to obtain the reflection times of the top and bottom interfaces of the target layer; based on this, combined with along-layer velocity coherence analysis, obtain the velocities of the overlying strata of the target layer and the velocity of the target layer, as shown in Figure 4 and Figure 5 respectively;

[0081] The above is only a preferred embodiment of the present invention and is not a limitation of the present invention in other forms. Any person skilled in the art may use the disclosed technical content to make changes or modifications into equivalent embodiments with equivalent changes and apply them to other fields. However, as long as it does not depart from the technical solution content of the present invention, any simple modification, equivalent change, and modification made to the above embodiments based on the technical essence of the present invention still fall within the protection scope of the technical solution of the present invention.

Claims

1. A method for identifying marine gas hydrates without typical seismic reflection characteristics, characterized in that, It includes the following steps: Step A: Optimize the pre-stack CMP gather; Step B: Horizon picking: On the stacked section, pick up the reflection times of the seabed and the top and bottom interfaces of the target horizons according to the seismic reflection characteristics; Step C: Along-layer coherent velocity analysis: On the optimized pre-stack CMP gather, with the seabed and horizon information picked up in Step B as constraints, obtain the velocity of the target horizon based on the coherent velocity inversion algorithm. Specifically: Step C1: Calculate the interval velocity using the seismic stacking velocity according to the Dix formula; ； In the formula, is the interval velocity between the -th and the -th interfaces, is the stacking velocity of the -th interface, is the zero-offset travel time of the -th interface; Step C2: For a specific horizon, set a series of trial velocity values with a certain velocity increment based on the initial velocity; Step C3: For a given trial velocity value, combine the reflection times of the top and bottom interfaces of the target horizon, calculate the depth of the interface, and calculate the travel time of the reflected wave in the CMP gather through ray tracing; Step C4: Calculate the coherence coefficient of the CMP gather along the calculated travel time curve; ； Where S l,m represents the coherence coefficient of the th CMP gather, represents the time window width, represents the number of seismic traces in the gather, represents the amplitude, is the offset of the th trace, time , (V n ) m represents the current test velocity; Step C5: Calculate the coherence coefficient corresponding to each trial velocity; the trial velocity corresponding to the maximum coherence coefficient is the interval velocity at the position of this CMP gather in the target horizon; Step C6: Process multiple CMP gathers for velocity picking according to Steps C1 - C5, obtain the interval velocities at the positions of multiple CMP gathers, interpolate them to obtain the interval velocities at all CMP points, and then perform smoothing processing to establish an initial model of the lateral variation of the interval velocity; Step C7: Calculate the coherence function of all CMP gathers according to the current initial velocity model; ； In the formula, is the coherence function, represents the number of CMP gathers; Step C8: Continuously iteratively adjust the velocity model. When the error between the coherence function calculated in this iteration and the coherence function calculated in the previous iteration is less than the error threshold, the velocity in this iteration is the result of the coherent velocity inversion, and thus obtain the velocity of the target horizon; Step D: Velocity feature extraction: Combine Step C to calculate the velocity difference between the target horizon and the overlying strata. Represent all negative values as velocity anomalies and positive values as normal velocities, and complete the calculation of the velocity differences for all survey lines; Step E: Velocity anomaly interpretation: Combine Step D to obtain the plane distribution map of the velocity differences of the target horizon in the work area. The velocity anomaly area corresponds to the hydrate distribution area.

2. The method for identifying marine hydrates without typical seismic reflection characteristics according to claim 1, wherein: The specific implementation of Step B is as follows: (1) On the stacked section, first automatically pick up the seabed reflection time according to the energy ratio method. For the case where the direct wave and the seabed reflection energy overlap, or the seabed isochron is deformed due to noise, combine the energy ratio method and the correlation method to determine the seabed reflection time; (2) Determine the top and bottom interfaces of the target horizon based on whether it is parallel to the seabed reflection isochron, and whether there is partial amplitude enhancement or weakening of the isochron, or the phenomenon of the isochron bulging upward.

3. The method for identifying marine hydrates without typical seismic reflection characteristics according to claim 1, wherein: The specific implementation of velocity feature extraction in Step D is as follows: If there is a high-speed hydrate thin layer immediately below the seabed, the seabed velocity shows a high-speed anomaly at the corresponding position, and the lateral distribution range of the high-speed anomaly indicates the lateral distribution range of the high-speed hydrate layer; for the target horizon below the seabed, calculate the velocity difference between the target horizon and the overlying strata. Represent all negative values as velocity anomalies and positive values as normal velocities.

4. The method for identifying marine hydrates without typical seismic reflection characteristics according to claim 1, wherein: The optimization operations in step A include noise suppression, flattening of common shot gathers, and excision of far offset stretch.

Images