A method for identifying the overlapping area between river channel sand bodies

By improving the time window pickup and data zero-complement method, combined with Ito conditions and fast Fourier transform, the problem of noise and time window selection is solved, the accurate identification of the superposition area of ​​the river phase sand body is achieved, and the resolution and calculation stability of the phase spectrum are improved.

CN115826055BActive Publication Date: 2025-07-04SOUTHWEST PETROLEUM UNIV
View PDF 0 Cites 0 Cited by

Patent Information

Application Number
CN202211715080.X
Authority / Receiving Office
CN · China
Patent Type
Patents(China)
Current Assignee / Owner
Filing Date
2022-12-28
Publication Date
2025-07-04
Estimated Expiration
2042-12-28

AI Technical Summary

Technical Problem

When identifying the overlapping area of ​​river phase sand bodies, the prior art is greatly affected by noise, time window selection, discrete sampling accuracy and seismic signal starting point, resulting in low phase spectrum resolution and disordered calculation results, making it difficult to accurately identify the overlapping area of ​​sand bodies.

Method used

The time window pickup method is improved. By using the first positive and negative transformation point outside the 'peak and valley extreme value' as the time window pickup point, and the seismic data is zeroed to 1001 points, combined with Ito conditions and fast Fourier transform, the unfolded phase spectrum is obtained, and integral calculation is performed to identify the superimposed area.

Benefits of technology

It improves the accuracy and stability of phase calculation, reduces noise interference, enhances the resolution of the phase spectrum, and realizes clear and intuitive recognition of the overlapping areas of river phase sand bodies.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN115826055B_ABST
    Figure CN115826055B_ABST
Patent Text Reader

Abstract

The present invention discloses a method for identifying the superimposed area between river channel sand bodies, which includes: inputting seismic signal data and performing time window picking on the seismic data of the target horizon; calculating the Nyquist cut-off frequency Nf of the phase spectrum according to the sampling rate of the seismic data, and zero-padding each trace of data; performing fast Fourier transform processing on the zero-padded seismic data to obtain the principal value phase spectrum of the seismic data; performing zero-padding processing on the phase data at the zero frequency of the obtained principal value phase spectrum; obtaining the expanded phase spectrum based on the Ito condition for the processed principal value phase spectrum; performing integral calculation on the expanded phase spectrum to obtain the integral expanded phase spectrum, and analyzing the integral expanded phase spectrum to identify the superimposed area of the river sand bodies. The advantages of the present invention are: improving the anti-noise interference ability, the accuracy of the phase calculation result, the resolution of the phase spectrum, facilitating the comparison and identification of the expanded phase, improving the stability of the result, and making the identification clearer and more intuitive.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] The present invention relates to the technical field of oil and gas exploitation, and particularly relates to a method for identifying the superimposed area between channel sand bodies based on an improved integral expansion phase spectrum attribute. Background Technique

[0002] Fluvial sand bodies are important oil and gas reservoirs in continental basins in China. Their sedimentary characteristics mainly include: uneven thickness vertically, with sand and mud interbedded; variable continuity horizontally, and rapid lithological changes. This is due to the frequent channel avulsion and migration, combined with factors such as diagenetic transformation and tectonic movement, resulting in the easy development of different types of lateral and stacked sand bodies in fluvial sand bodies, forming thin interbedded superimposed areas of mudstone-sandstone-mudstone. The spatial distribution of such thin interbedded sand bodies and the identification of these superimposed areas are one of the main research objectives of seismic exploration in oil and gas fields, and are of great significance for the study of remaining oil distribution in an old oilfield, as well as the deployment and adjustment of development plans.

[0003] Currently, in the identification of the superimposed area of fluvial sand bodies, seismic attribute slices and stratigraphic slices based on seismic sedimentology are usually used to identify sand body superimposition, or a combination of multiple attributes such as instantaneous attributes and along-layer average attributes is used to identify the superimposed area. For the exploitation of oil and gas reservoirs in fluvial sand body reservoirs, most are faced with thin reservoirs below the seismic wave tuning thickness, and these thin reservoirs usually show a waveform of one cycle in seismic records. For such a waveform of one cycle, the changes in amplitude and frequency are often used to reflect the changes in thin layer thickness, but it is difficult to directly identify the seismic waveform changes caused by the superimposition of two channel sand bodies. In the application of seismic wave phase, it is mainly based on the principal value phase after Fourier transform (within a 2π period) and the instantaneous phase and its expansion phase based on Hilbert transform. Currently, it is mainly applied to formation thickness estimation, seismic horizon tracking, and the identification of faults, stratigraphic unconformities, etc., and there is little application of phase to the identification of the superimposed area between channel sand bodies.

[0004] Liu Yang (2019 [1] ) and Wang Banni (2020 [2] , 2022 [3] ), on the basis of summarizing the previous research on the prediction method of fluvial discontinuity boundaries and the application of phase spectrum, considering that geological structure changes are more likely to affect the phase characteristics of seismic reflection waves, proposed a method using integral expansion phase spectrum to identify the discontinuity boundaries of fluvial sand bodies. Through theoretical model and actual data tests, it demonstrated certain potential of integral expansion phase spectrum in the identification of fluvial superimposed sand bodies;

[0005] However, the above-mentioned existing technologies are greatly affected by the selection of signal time window and noise level, and have the following specific defects:

[0006] First, since seismic signal data always contains a large amount of various types of noise, the noise will affect the calculation of the phase when extracting the phase attributes of the seismic waves, and the noise will also affect the picking effect of the time window of the seismic signal, thereby affecting the effect of identifying the sand-shale superimposed area by the integral phase spectrum;

[0007] Second, due to the reason of discrete sampling accuracy, it is usually difficult to find a zero point outside the "peak-valley extreme value". The found "zero point" may be any point around the waveform zero point, and the position of the "zero point" of each trace may be different, resulting in the chaotic seismic signals picked by the time window, thereby affecting the subsequent calculation;

[0008] In addition, since the seismic signals picked by the time window are short, and the lengths of the signals picked by different traces vary due to the change of sand body thickness, it will also affect the resolution of the phase spectrum;

[0009] Fourth, since the starting point of the extracted seismic signal is not the true zero point of the waveform, but positive or negative at the starting point of the time window, the starting point (zero frequency point) of the calculated principal value phase spectrum will also become chaotic, thus affecting the calculation of the subsequent unfolded phase spectrum.

[0010] References:

[0011] [1] Liu Yang. Application of Unfolded Phase Spectrum in the Discontinuity Boundary of Superimposed Sand Bodies [D]. Southwest Petroleum University, 2019;

[0012] [2] Wang Banni, Yin Cheng, Liu Yang, et al. Research on the Identification Method of the Discontinuity Boundary between Superimposed Sand Bodies Based on Unfolded Phase Spectrum [A]. 2020 Annual Meeting of Chinese Geoscience Union [C], 2020;

[0013] [3] Wang Banni, Yin Cheng, Liu Yang, et al. Identification of Sand Body Discontinuity Boundary Based on Unfolded Phase Spectrum [J]. China Petroleum and Chemical Standard and Quality, 2022, 42(1): 165-167.

[0014] Definitions and Explanations of Key Terms:

[0015] Sand body superimposition: Due to the change of sedimentary environment, the channel sand bodies with semi-elliptical forms often migrate and change their courses, so that the channel sand bodies in different periods may overlap partially in the vertical and horizontal directions, which is called sand body superimposition.

[0016] Sand body superimposition area: Two phases of channel sand bodies deposited in the same area at different times can either have adjacent boundaries or overlap by cutting a certain range. Due to changes in factors such as accommodation space and sediment source rate in the vertical direction, the top interfaces of the two channel sand bodies may not be at the same horizontal plane, resulting in a certain elevation difference between the two phases of channels. Therefore, the superimposition of the two phases of channel sand bodies may present as an area with different widths and thicknesses, which is called the sand body superimposition area.

[0017] Discontinuity boundary: For two phases of channel sand bodies in the sand body superimposition area, due to differences in deposition time and environment, their physical properties often vary. Additionally, in the superimposition area, if the two phases of channels are not in complete close contact, some dense mudstones often fill the middle. Therefore, this local physical property heterogeneity change in the sand body reservoir is called the discontinuity boundary of the reservoir.

[0018] Principal value phase spectrum: According to Fourier transform theory, for the spectrum of a signal, the phase magnitudes (i.e., the positions of vibrations) of multiple harmonics of this signal from low frequency to high frequency can be directly obtained using the real and imaginary parts of the Fourier transform through the arctangent or arcsine function. Due to the periodicity of the sine function, this phase spectrum generally shows values between 0 and 2π or -π and +π, so it is called the principal value phase spectrum of this signal.

[0019] Unwrapped phase spectrum: The position of a signal's vibration (i.e., the phase) can increase from 0 to π / 2, π, 3π / 2, 2π, and continue to increase to infinity. Therefore, to overcome the problem of periodic folding in obtaining the phase spectrum by Fourier transform, it is necessary to unwrap the conventional principal value phase spectrum. Thus, the phase spectrum obtained by using the phase unwrapping method based on the principal value phase is called the unwrapped phase spectrum. From a mathematical perspective, phase unwrapping is an ill - posed problem, and currently, various algorithms approximate the true solution to a certain extent.

[0020] Integrated unwrapped phase spectrum: To highlight the differences in the unwrapped phase spectra between different seismic traces, the unwrapped phase spectra of different frequencies obtained from the same seismic trace are integrated and summed to obtain an accumulated phase value, and the accumulated phase values obtained from all seismic traces are arranged together, which is called the integrated unwrapped phase spectrum. Summary of the Invention

[0021] Aiming at the defects of the prior art, the present invention provides a method for identifying the superimposition area between channel sand bodies. The window - picking method is improved to reduce the influence of noise. At the same time, the methods for obtaining the principal value phase and the unwrapped phase are improved, making the identification effect clearer and more intuitive.

[0022] To achieve the above - mentioned invention purposes, the technical solutions adopted by the present invention are as follows:

[0023] A method for identifying the superimposed area between channel sand bodies, the specific steps are as follows:

[0024] Step 1: Input seismic signal data and perform time window picking on the seismic data of the target horizon. In order to include a complete composite wave of thin interbeds to the greatest extent, discard the interference of redundant reflection waves from other horizons, and eliminate the ineffective identification and chaotic influence of the zero points outside the "peak-valley extreme value" of the "one-period time window", the picking method of the "one-period time window" is improved. The first positive-negative transformation point outside the "peak-valley extreme value" is used as the point for time window picking, and in order to weaken the influence of noise, the distance for searching for "zero points" by extending the "peak-valley extreme value" is restricted according to the seismic data;

[0025] Step 2: Since the sampling rate of the seismic data is 1 ms, calculate the Nyquist cut-off frequency Nf of the phase spectrum to be 1000 Hz. Therefore, each trace of data is padded with zeros to 1001 points, so that the phase spectrum data in the frequency domain has corresponding phase change data at each 1 Hz frequency in the range from 0 to 1000 Hz.

[0026] Step 3: Perform fast Fourier transform processing on the seismic data padded with zeros in Step 2 to obtain the principal value phase spectrum of the seismic data;

[0027] Step 4: Pad the phase data at the zero frequency of the principal value phase spectrum obtained in Step 3 with zeros to improve the stability of the calculation result; the zero-padding process is to replace all the chaotic phase data at the zero frequency in the principal value phase spectrum with zero-phase data.

[0028] Step 5: Based on the Ito condition, obtain the expanded phase spectrum for the principal value phase spectrum processed in Step 4

[0029] Step 6: Perform integral calculation on the obtained expanded phase spectrum to further highlight the differences in the expanded phase spectra between different seismic traces, which is called the integral expanded phase spectrum,

[0030] Step 7: Analyze the integral expanded phase spectrum to identify the superimposed area of river sand bodies. In the integral expanded phase spectrum obtained after processing, the non-superimposed area of thin interbeds shows a lower background value, and the superimposed area shows an abnormally high value. Identify the sand body superimposed area according to the abnormal mutation position of the graph.

[0031] Furthermore, in Step 2, the calculation formula of Nf is as follows:

[0032]

[0033] where, △t is the sampling interval of the seismic signal, such as 1 ms mentioned above.

[0034] Furthermore, the sub-steps of Step 5 are as follows:

[0035] S1: First, determine the position of the unwrapping point of the unwrapped phase. This can be done by judging the phase difference Δθ(ω i , Δω) between adjacent frequency sampling points to determine the position of the unwrapping point, that is

[0036] Δθ(ω i , Δω) = θ p (ω i + Δω) - θ p (ω i )

[0037] where Δω = ω i+1 - ω i , and i is the serial number of the frequency sampling point.

[0038] S2: According to the Itô condition, the phase difference between adjacent frequency sampling points should satisfy:

[0039] - π < Δθ(ω i , Δω) ≤ π

[0040] If the phase difference does not satisfy the Itô condition, it is considered a discontinuous point at this time, and thus the phase spectrum is unwrapped. To eliminate the discontinuity phenomenon, we usually add an appropriate 2kπ to the conventional phase spectrum. That is

[0041]

[0042] where k is the number of times of unwrapping the phase, starting from 1 and increasing in sequence. Among them, is the unwrapped phase spectrum, and ω is the frequency.

[0043] Furthermore, in step six, the formula for obtaining the integral unwrapped phase spectrum is as follows:

[0044]

[0045] where θ H is the integral unwrapped phase spectrum, is the unwrapped phase spectrum, and ω is the frequency.

[0046] Compared with the prior art, the advantages of the present invention are as follows:

[0047] First, improve the method of time window extraction. It not only unifies the starting position of the signals picked up by the time window, eliminates the influence of the disordered zeros between channels, but also improves the anti-interference ability against noise, effectively picks up the seismic records in the sand body stacking area, and thus improves the accuracy of the phase calculation result;

[0048] Secondly, improve the method for obtaining the principal value phase. Before calculating the principal value phase, zero-padding the extracted short-time window seismic data to 1001 points improves the resolution of the phase spectrum and unifies the data lengths of all channels, facilitating the comparison and identification of the unwrapped phases.

[0049] Thirdly, improve the method for obtaining the unwrapped phase. Before calculating the unwrapped phase, zero-padding the phase at the zero frequency of the principal value phase spectrum improves the stability of the calculation result of the unwrapped phase spectrum, making the identification of the stacked area of fluvial sand bodies using the integral unwrapped phase spectrum clearer and more intuitive. Description of the Drawings

[0050] Figure 1 is the flowchart of the method for identifying the stacked area between channel sand bodies in the embodiment of the present invention;

[0051] Figure 2 is the schematic diagram of time window selection in the embodiment of the present invention;

[0052] Figure 3 is the diagram of four types of wedge-shaped sand body models in the embodiment of the present invention; among them, the light color is sandstone and the dark color is mudstone;

[0053] Figure 4 is the unwrapped phase spectrum of four time windows in the embodiment of the present invention; (a) Unwrapped phase spectrum of AC time window, (b) Unwrapped phase spectrum of AD time window, (c) Unwrapped phase spectrum of BC time window, (d) Unwrapped phase spectrum of BD time window;

[0054] Figure 5 is the integral phase spectrum of four time windows in the embodiment of the present invention; (a) Unwrapped phase spectrum of AC time window, (b) Unwrapped phase spectrum of AD time window, (c) Unwrapped phase spectrum of BC time window, (d) Unwrapped phase spectrum of BD time window;

[0055] Figure 6 is the five types of wedge-shaped sand body models in the embodiment of the present invention; among them, the light color is sandstone and the dark color is mudstone;

[0056] Figure 7 is the effect diagram of the time window picking trajectory before and after improvement in the embodiment of the present invention (the thick solid line is the start point of the time window; the thin solid line is the end point of the time window); (a) Time window picking trajectory without noise (before improvement), (b) Time window picking trajectory without noise (after improvement), (c) Time window picking trajectory with 5% noise (before improvement), (d) Time window picking trajectory with 5% noise (after improvement), (e) Time window picking trajectory with 10% noise (before improvement), (f) Time window picking trajectory with 10% noise (after improvement), (g) Time window picking trajectory with 20% noise (before improvement), (h) Time window picking trajectory with 20% noise (after improvement);

[0057] Figure 8It is the integral expansion phase contrast with different numbers of zero-padding in the AD time window of the embodiment of the present invention. (a) Without zero-padding, (b) Zero-padding to 101 points, (c) Zero-padding to 401 points, (d) Zero-padding to 1001 points, (e) Zero-padding to 1201 points, (f) Zero-padding to 2001 points;

[0058] Figure 9 It is a model of seven types of wedge-shaped sand bodies in the embodiment of the present invention; among them, the light color is sandstone and the dark color is mudstone;

[0059] Figure 10 It is a comparison diagram of the expanded phase spectra before and after zero-padding at zero frequency in the embodiment of the present invention. (a) Expanded phase spectrum before zero-padding, (b) Expanded phase spectrum after zero-padding;

[0060] Figure 11 It is a comparison diagram of the integral expanded phase spectra before and after zero-padding at zero frequency in the embodiment of the present invention; (a) Integral expanded phase spectrum before zero-padding, (b) Integral expanded phase spectrum after zero-padding. Detailed implementation manners

[0061] To make the objectives, technical solutions and advantages of the present invention clearer and more understandable, the following further elaborates on the present invention in detail according to the attached drawings and by listing embodiments.

[0062] As Figure 1 shown, a method for identifying the superimposed area between channel sand bodies is as follows:

[0063] Step 1: Input seismic signal data and pick up the time window of the seismic data of the target horizon. In order to include a complete composite wave of thin interbeds to the greatest extent, discard the interference of redundant reflection waves of other horizons, and eliminate the ineffective identification and chaotic influence of the zero points outside the "peak-valley extreme value" of "one-period time window", the picking method of "one-period time window" is improved. As Figure 2 shown, the first positive-negative transformation points outside the "peak-valley extreme value" (i.e., point A and point D) are used as the points for picking up the time window, and in order to weaken the influence of noise, the distance for searching for "zero points" by extending the "peak-valley extreme value" is restricted according to the seismic data;

[0064] Step 2: Zero-pad the tail end of the picked finite-length seismic data to 1001 points. After zero-padding, the corresponding phase data changes at each 1 Hz frequency are realized; the lengths of zero-padding for seismic traces with different data lengths are different, and finally all trace data are zero-padded to the same value.

[0065] According to the sampling rate of the seismic data being 1 ms, the Nyquist cut-off frequency (Nf) of the phase spectrum is calculated to be 1000 Hz. Therefore, each trace data is zero-padded to 1001 points, so that the phase spectrum data in the frequency domain has corresponding phase change data at each 1 Hz frequency in the range from 0 to 1000 Hz. The calculation formula is as follows:

[0066]

[0067] Step 3: Perform fast Fourier transform processing on the seismic data after zero-padding in Step 2 to obtain the principal value phase spectrum of the seismic data;

[0068] Step 4: Perform zero-padding on the phase data at the zero frequency of the principal value phase spectrum obtained in Step 3 to improve the stability of the calculation results;

[0069] Zero-padding processing: All the chaotic phase data at the zero frequency in the principal value phase spectrum are replaced with zero phase data.

[0070] Step 5: Obtain the expanded phase spectrum based on the Itô condition for the principal value phase spectrum processed in Step 4

[0071] First, determine the position of the expansion point of the expanded phase. This can be determined by judging the phase difference Δθ(ω i ,Δω) between adjacent frequency sampling points to determine the position of the expansion point, that is

[0072] Δθ(ω i ,Δω) = θ p (ω i +Δω) - θ p (ω i ) (1)

[0073] where Δω = ω i+1 -ω i , and i is the serial number of the frequency sampling point.

[0074] According to the Itô condition, the phase difference between adjacent frequency sampling points should satisfy:

[0075] -π < Δθ(ω i ,Δω) ≤ π (2)

[0076] If the phase difference does not satisfy the Itô condition, it is considered a discontinuous point at this time, and thus the phase spectrum is expanded. To eliminate the discontinuity phenomenon, we usually add an appropriate 2kπ to the conventional phase spectrum. That is

[0077]

[0078] where k is the number of times of the expanded phase, starting from 1 and increasing in sequence. Among them, is the expanded phase spectrum, and ω is the frequency.

[0079] Step 6: Perform integral calculation on the obtained expanded phase spectrum to obtain the integral expanded phase spectrum;

[0080] Integrating the unwrapped phase spectra of different frequencies to further highlight the differences in the unwrapped phase spectra between different seismic traces is called the integrated unwrapped phase spectrum, and the formula is as follows:

[0081]

[0082] where θ H is the integrated unwrapped phase spectrum, is the unwrapped phase spectrum, and ω is the frequency.

[0083] Step 7: Analyze the integrated unwrapped phase spectrum to identify the superimposed areas of river sand bodies.

[0084] The integrated unwrapped phase spectrum obtained after the above processing is as shown in Figure 11 (b). The non-superimposed areas of thin interbeds show lower background values, and the superimposed areas show abnormally high values. The sand body superimposed areas can be identified according to the abnormal mutation positions of the graph.

[0085] Improving the time window picking method can more effectively extract seismic records, making the obtained unwrapped phase spectrum clearer. In addition, the improved time window picking method has stronger anti-noise interference ability, reducing the interference of noise on the identification of superimposed areas and increasing the allowable noise content in the sandstone-shale interbeds that can be identified.

[0086] Since the zero points outside the "peak-valley extreme values" of the "one-period time window" cannot be effectively identified and are chaotic, as shown in Figure 2 , the first "zero point" outside the "peak-valley extreme values" can be unified into two points near the true zero point, and four time window picking schemes can be formed, namely AC, AD, BC, and BD.

[0087] For Figure 3 the wedge-shaped sand body model, the impacts on subsequent phase calculations of the time windows picked according to the four different schemes are as shown in Figure 4 and Figure 5 :

[0088] Among the four schemes, picking the AD time window for subsequent calculations results in a higher accuracy rate in identifying the superimposed areas in the unwrapped phase spectrum and the integrated phase spectrum. Thus, it can be determined that the "period time window" determined by AD is the most effective time window picking method for identifying the superimposed areas.

[0089] To weaken the interference of noise, in the present invention, certain restrictions are made based on seismic data when determining the four points A, B, C, and D outside the "peak-valley extreme values". Taking the model of Figure 6 as an example, the time window picking effect is as shown in Figure 7 :

[0090] The improved window picking method has stronger anti-interference ability against noise, reduces the interference of noise on stacked identification, and improves the anti-noise ability when identifying sandstone-shale interbeds.

[0091] Zero-padding processing is performed on the seismic records extracted by the window picking method. Through multiple experiments, it is confirmed that zero-padding to 1001 points enables a corresponding phase change for each 1 Hz frequency change, which can make the resolution of the calculated phase data the highest and the effect the best. Similarly, taking Figure 4 four types of wedge models as an example, when the window data is zero-padded to different lengths, the comparison of the integral expanded phase spectra is as Figure 8 shown:

[0092] Before zero-padding to 1001 points, the more the number of zero-padding points, the better the performance in terms of details. After zero-padding exceeds 1001 points, the detail performance is worse than that at 1001 points;

[0093] Phase zero-padding processing is performed at zero frequency on the obtained principal value phase spectrum (assigning zero to the first sampling point of the principal value phase of all traces), which can improve the stability of the calculation result of the expanded phase spectrum and make the expanded phase data of all seismic traces more regular.

[0094] Taking Figure 9 seven types of wedge models as an example, through Figure 10 and Figure 11 it can be obtained that after zero-phase zero-padding processing on the principal value phase spectrum, the stability of the calculated expanded phase spectrum can be improved, making the expanded phase data of all seismic traces more regular, and the integral expanded phase spectrum is also easier to identify the stacked area.

[0095] Those of ordinary skill in the art will realize that the embodiments described herein are for helping readers understand the implementation methods of the present invention, and should be understood that the protection scope of the present invention is not limited to such specific statements and embodiments. Those of ordinary skill in the art can make various other specific deformations and combinations that do not depart from the essence of the present invention based on the technical revelations disclosed in the present invention, and these deformations and combinations are still within the protection scope of the present invention.

Claims

1. A method for identifying the superimposed area between river channel sand bodies, characterized in that, The specific steps are as follows: Step 1: Input seismic signal data, and perform time window picking on the seismic data of the target horizon; improve the "one-cycle time window" picking method; take the first positive-negative transformation point outside the "peak-valley extreme value" as the time window picking point, and in order to weaken the influence of noise, make a limit on the distance of searching for the "zero point" by extending the "peak-valley extreme value" according to the seismic data; Step 2: Since the sampling rate of the seismic data is 1 ms, calculate the Nyquist cut-off frequency Nf of the phase spectrum to be 1000 Hz. Therefore, zero-padding each trace of data to 1001 points makes the phase spectrum data in the frequency domain have corresponding phase change data at each 1 Hz frequency within the range from 0 to 1000 Hz; Step 3: Perform fast Fourier transform processing on the zero-padded seismic data in Step 2 to obtain the principal value phase spectrum of the seismic data; Step 4: Perform zero-padding processing on the phase data at the zero frequency of the principal value phase spectrum obtained in Step 3 to improve the stability of the calculation result; The zero-padding processing is to replace all the chaotic phase data at the zero frequency in the principal value phase spectrum with zero-phase data; Step 5: Obtain the expanded phase spectrum based on the Itô condition for the principal value phase spectrum processed in Step 4; Step 6: Perform integral calculation on the obtained expanded phase spectrum to further highlight the differences in the expanded phase spectra between different seismic traces, which is called the integral expanded phase spectrum; Step 7: Analyze the integral expanded phase spectrum to identify the superimposed areas of fluvial sand bodies; in the integral expanded phase spectrum obtained after processing, the non-superimposed areas of thin interbeds show lower background values, and the superimposed areas show abnormally high values. Identify the sand body superimposed areas according to the abnormal mutation positions of the graph.

2. The method for identifying the superimposed area between river channel sand bodies according to claim 1, wherein: In Step 2, the calculation formula of Nf is as follows: where, △t is the sampling interval of the seismic signal, such as 1 ms mentioned above.

3. A method for identifying the superimposed area between river channel sand bodies according to claim 1, characterized in that: The sub-steps of Step 5 are as follows: S1: First, determine the position of the unwrapping point of the unwrapped phase. This can be done by judging the phase difference Δθ(ω i ,Δω) of adjacent frequency sampling points to determine the position of the unwrapping point, that is Δθ(ω i ,Δω) = θ p (ω i +Δω) - θ p (ω i ) where Δω = ω i+1 - ω i , and i is the serial number of the frequency sampling point; S2: According to the Itô condition, the phase difference between adjacent frequency sampling points should satisfy: -π < Δθ(ω i , Δω) ≤ π When the phase difference does not satisfy the Itô condition, it is considered as a discontinuous point at this time, and thus the phase spectrum is expanded; in order to eliminate the discontinuous phenomenon, we usually add an appropriate 2kπ to the conventional phase spectrum; that is where k is the number of times of the unwrapped phase, increasing sequentially from 1; where is the unwrapped phase spectrum, and ω is the frequency.

4. A method for identifying the superimposed area between channel sand bodies according to claim 1, characterized in that: In Step 6, the formula for obtaining the integral expanded phase spectrum is as follows: Among them, θ H is the integral expansion phase spectrum, is the expansion phase spectrum, and ω is the frequency.