Seismic weak signal detection method, system and device based on instantaneous phase analysis
The seismic weak signal detection method based on instantaneous phase analysis solves the problem of difficulty in identifying weak reflected signals from thin reservoirs or target layers, enabling more accurate reservoir prediction and stratigraphic interpretation.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- PETROCHINA CO LTD
- Filing Date
- 2024-11-25
- Publication Date
- 2026-05-26
Smart Images
Figure CN122085337A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of reservoir prediction technology in geophysical exploration of oil and gas, and particularly to a method, system, and equipment for detecting weak seismic signals based on instantaneous phase analysis. Background Technology
[0002] In recent years, with the increasing sophistication of oil and gas exploration, the focus has shifted towards thin reservoir exploration, often requiring reservoir thicknesses as thin as 10 meters or even lower. Therefore, detecting weak signals carrying seismic wave reflection information from thin reservoirs or target layers has become a major challenge in reservoir prediction within oil and gas geophysical exploration. The formation mechanism of weak signals carrying effective information in seismic reflection waves is that the rock physical properties and wave impedance of thin reservoirs or target layers are relatively similar to those of the overlying strata. Consequently, the energy or amplitude of reflected seismic waves from such thin reservoirs or target layers is weak, making accurate identification difficult and thus a significant obstacle to accurate reservoir prediction.
[0003] Currently, there are three main types of technologies in the field of weak seismic signals: the first type is weak signal extraction methods based on compressed sensing theory; the second type is weak signal energy enhancement methods based on spectral enhancement, such as those based on generalized S-transform, wavelet transform, maximum entropy spectrum method, and curvelet transform; and the third type is methods based on instantaneous attribute analysis, such as those based on Hilbert transform and empirical mode decomposition.
[0004] However, careful analysis reveals that the aforementioned methods are generally ineffective in identifying and characterizing weak seismic signals. For example, Wei Xueqiang (2017) used S-transform and wavelet transform, but this method is designed for time-frequency analysis and requires signals to have a certain degree and scale of response differences, rather than being specifically designed for weak seismic signals. Qu Zhongdang (2015) used S-transform to extract the difference between weak signals and noise in the frequency domain, but this suppressed mixing interference within the effective frequency band. Jin Dan (2016) used curvelet transform to identify the difference in correlation between weak signals and noise, thereby removing noise and highlighting the effective signal, but this method primarily aims to identify noise. Similarly, Zhou Yang (2018) used compressed sensing, and Le Youxi (2024) used K-means singular value decomposition to identify weak signals, but their methods focus on noise identification and removal. Therefore, the above methods are less effective when applied to the detection and identification of weak seismic signals and cannot provide accurate technical support for seismic horizon interpretation or reservoir prediction. Summary of the Invention
[0005] The purpose of this invention is to overcome the aforementioned problems in the prior art and provide a seismic weak signal detection method, system, and device based on instantaneous phase analysis. This detection method first performs phase shifting on the seismic data to avoid the influence of Hilbert transform on the phase; then, it decomposes the seismic data into a series of common-frequency data volumes based on time-frequency decomposition; next, it solves for the second-order instantaneous phase based on Hilbert transform to obtain a series of three-dimensional weak signal indicator data volumes; finally, it selects the data volume with the best matching relationship with the weak signal of the target layer in the study area as the optimal three-dimensional weak signal indicator data volume, thereby establishing a correspondence between weak signals and characteristic data volumes, achieving the effect of effectively identifying and indicating weak signal reflections, and providing technical support and basis for seismic horizon interpretation or reservoir prediction.
[0006] To achieve the above objectives, the technical solution adopted by the present invention is as follows:
[0007] In a first aspect, the present invention provides a method for detecting weak seismic signals based on instantaneous phase analysis, which includes the following steps:
[0008] Step S1: Input the 3D seismic data volume of the study area, as well as the top interface depth data and time-depth relationship data of the target layer after drilling, and denote the dominant frequency of the seismic signal of the target layer in the 3D seismic data volume as Fre_main; wherein, the 3D seismic data volume includes the x-axis, y-axis and z-axis, the x-axis represents the direction of the survey line, and the number of lines of the survey line is denoted as N_line; the y-axis represents the direction of the seismic trace, and the total number of seismic traces of a single survey line is denoted as N_trace; the z-axis represents time, and the total number of time sampling points of a single seismic trace is denoted as N_time;
[0009] Step S2: Set the current seismic trace line number as id_line and the seismic trace number as id_trace; extract the amplitude data of the current seismic trace from the three-dimensional seismic data volume and record it as a one-dimensional array x(n), the length of which is equal to the total number of time sampling points N_time of the current seismic trace;
[0010] Step S3: Perform a -90° phase shift on the one-dimensional array x(n) to obtain the phase-shifted array y(n);
[0011] Step S4: Perform time-frequency continuous wavelet transform on the phase shift array y(n) to obtain multiple common frequency amplitude arrays YAi(n) in the mid-to-high frequency band;
[0012] Step S5: Perform Hilbert transform on each common frequency amplitude array YAi(n) to obtain the instantaneous phase array YA_Phasei(n) corresponding to each common frequency amplitude array YAi(n);
[0013] Step S6: For each instantaneous phase array YA_Phasei(n), calculate the second derivative Gi(n) with respect to the time variable n to obtain the three-dimensional data volume of the weak signal indication of the current seismic trace co-frequency.
[0014] Step S7: Change the seismic trace and repeat steps S2-S6 until all seismic traces in the study area are calculated, and obtain multiple weak signal indicators of the same frequency in the mid-to-high frequency band for three-dimensional data volume.
[0015] Step S8: Match all three-dimensional data volumes of weak signals with the same frequency with the top interface depth data and time-depth relationship data input in step S1. Select the three-dimensional data volume of weak signals with the same frequency that has the highest time-depth relationship with the weak reflection interface of the target section of the completed well in the study area. Record it as the optimal weak signal indicator data volume of the study area and complete the detection.
[0016] In step S3, the process of performing a -90° phase shift on the one-dimensional array x(n) is as follows:
[0017] Step S3.1: Perform a discrete Fourier transform on the one-dimensional array x(n) to obtain the spectrum array X(k), where K∈[0,N_time-1], and the part of k∈[(N_time-1) / 2,N_time-1] corresponds to the negative frequency;
[0018] Step S3.2: Perform frequency division calculation on the spectrum array X(k) to obtain the frequency-divided array Y(k);
[0019] Step S3.3: Perform an inverse discrete Fourier transform on the frequency-division processed array Y(k), take its real part, and obtain the phase-shifted array y(n).
[0020] In step S3.2, the formula for calculating the frequency division of the spectrum array X(k) is as follows:
[0021]
[0022] In step S4, the frequency range of the common frequency amplitude array YAi(n) is from 0.8 times the main frequency to 2 times the main frequency; where YA1(n) represents the amplitude array with a frequency of 0.8×Fre_main, YA2(n) represents the amplitude array with a frequency of 0.8×Fre_main+Δf, YA3(n) represents the amplitude array with a frequency of 0.8×Fre_main+2×Δf, ..., Δf takes the value 2 to 5.
[0023] In step S5, the process of performing a Hilbert transform on the common frequency amplitude array YAi(n) is as follows:
[0024] Step S5.1: Perform a discrete Fourier transform on the common frequency amplitude array YAi(n) to obtain the spectrum array YA_Fi(k), where K∈[0,N_time-1], and the part of k∈[(N_time-1) / 2,N_time-1] corresponds to the negative frequency;
[0025] Step S5.2: Perform frequency division calculations on the spectrum array YA_Fi(k) to obtain the frequency-divided array YT. i (k);
[0026] Step S5.3: Divide the frequency-division array YT i (k) Perform inverse discrete Fourier transform to obtain the analytic signal array YA_HTi(n). Record the real part data in the analytic signal array YA_HTi(n) as the real part array R_YA_HTi(n), and record the imaginary part data in the analytic signal array YA_HTi(n) as the imaginary part array I_YA_HTi(n).
[0027] Step S5.4: Calculate the instantaneous phase array YA_Phasei(n) based on the real part array R_YA_HTi(n) and the imaginary part array I_YA_HTi(n).
[0028] In step S5.2, the formula for calculating the frequency division of the spectrum array YA_Fi(k) is as follows:
[0029]
[0030] In step S5.4, the formula for calculating the instantaneous phase array YA_Phasei(n) is as follows:
[0031]
[0032] Secondly, the present invention provides a seismic weak signal detection system based on instantaneous phase analysis, comprising:
[0033] The data input module is used to input the three-dimensional seismic data volume of the study area, to input the top interface depth data and time-depth relationship data of the target layer after drilling, and to record the dominant frequency of the seismic signal of the target layer in the three-dimensional seismic data volume as Fre_main;
[0034] The setting and extraction module is used to set the current seismic trace's survey line number as id_line and the seismic trace number as id_trace; it is used to extract the amplitude data of the current seismic trace and record it as a one-dimensional array x(n), the length of which is equal to the total number of time sampling points N_time of the current seismic trace;
[0035] The phase shift processing module is used to perform a -90° phase shift processing on the one-dimensional array x(n) and obtain the phase shift array y(n);
[0036] The transformation module is used to perform time-frequency continuous wavelet transform on the phase shift array y(n) and obtain multiple common frequency amplitude arrays YAi(n) in the mid-to-high frequency band.
[0037] The Hilbert transform module is used to perform Hilbert transform on each common frequency amplitude array YAi(n) and obtain the instantaneous phase array YA_Phasei(n) corresponding to each common frequency amplitude array YAi(n);
[0038] The second derivative calculation module is used to calculate the second derivative Gi(n) with respect to the time variable n for each instantaneous phase array YA_Phasei(n), and obtain the three-dimensional data volume of weak signal indication of the current seismic trace co-frequency;
[0039] The integrated processing module is used to replace the seismic traces. It controls the setting and extraction module, phase shift processing module, transformation module, Hilbert transformation module and second derivative calculation module to repeatedly process the replaced seismic traces until all seismic traces in the study area are calculated and multiple weak signal indicators in the mid-to-high frequency band are obtained.
[0040] The optimal matching module is used to match all three-dimensional data volumes of weak signals with the same frequency with the top interface depth data and time-depth relationship data in the data input module. It is used to select the three-dimensional data volume of weak signals with the same frequency as the weak reflection interface time-depth relationship of the target section of the completed well in the study area with the highest consistency rate, thereby obtaining the optimal weak signal indicator data volume of the study area.
[0041] Thirdly, the present invention provides an electronic device comprising:
[0042] At least one processor; and,
[0043] A memory communicatively connected to the at least one processor; wherein,
[0044] The memory stores a computer program that can be executed by the at least one processor, which enables the at least one processor to perform the aforementioned method for detecting weak seismic signals based on instantaneous phase analysis.
[0045] Fourthly, the present invention provides a storage medium storing instructions that, when executed on a computer, cause the computer to perform the aforementioned method for detecting weak seismic signals based on instantaneous phase analysis.
[0046] The beneficial effects of using the present invention are as follows:
[0047] 1. In the detection method of the present invention, step S3 has the advantage of overcoming the inherent defects of the Hilbert transform, making the final data volume more accurate in responding to weak reflection signals. Step S4 has the advantage of selecting only the mid-to-high frequency band data from the frequency division processing results, relatively reducing the difficulty of data processing and accelerating the data processing speed. Step S6 has the advantage of further strengthening the indication ability of the instantaneous phase for weak signals by performing second-order differentiation on the instantaneous phase, which significantly improves and enhances the indication ability of the final three-dimensional data volume for weak signals. Step S6 has the advantage of comparing and optimizing multiple three-dimensional data volumes for weak signals from different frequency bands, thereby selecting the data volume most sensitive to weak signals, making the detection results basically consistent with the weak signals, making the detection results more reliable, and avoiding the common problem of multiple data volumes after frequency division processing without selection criteria.
[0048] In summary, this invention first performs phase shifting on seismic data to avoid the influence of Hilbert transform on the phase; then, it decomposes the seismic data into a series of common-frequency data volumes based on time-frequency decomposition; next, it solves for the second-order instantaneous phase based on Hilbert transform to obtain a series of three-dimensional data volumes indicating weak signals; finally, it selects the data volume with the best matching relationship with the weak signals of the target layer in the study area as the optimal three-dimensional data volume indicating weak signals, establishing a correspondence between weak signals and characteristic data volumes, thereby achieving the technical effect of effectively identifying and indicating weak signal reflections, providing technical support and basis for seismic horizon interpretation or reservoir prediction.
[0049] 2. This invention is based on seismic data reconstruction using instantaneous phase. Leveraging the sensitivity of instantaneous phase to geological structure, it establishes a correspondence between weak signals and characteristic data volumes through techniques such as phase shifting, frequency division, data optimization, and second derivative calculation. This effectively identifies three-dimensional data volumes indicating weak signals, improving the accuracy of seismic data reconstruction. Attached Figure Description
[0050] Figure 1 This is a flowchart of the present invention;
[0051] Figure 2 This is a system block diagram of the present invention;
[0052] Figure 3 This is the original seismic profile of two completed wells in a study area of Sichuan Province in Example 4;
[0053] Figure 4 For corresponding Figure 3 Instantaneous phase profile;
[0054] Figure 5This is the optimal weak signal indication profile obtained based on the present invention in Example 4. Detailed Implementation
[0055] Example 1
[0056] like Figure 1 As shown, this invention provides a method for detecting weak seismic signals based on instantaneous phase analysis, which includes the following steps:
[0057] Step S1: Input the 3D seismic data volume of the study area, as well as the top interface depth data and time-depth relationship data (in meters) of the target layer after drilling, and denote the dominant frequency of the seismic signal of the target layer in the 3D seismic data volume as Fre_main. The 3D seismic data volume includes x-axis, y-axis, and z-axis. The x-axis represents the direction of the seismic line, and the number of lines is denoted as N_line; the y-axis represents the direction of the seismic trace, and the total number of traces on a single seismic line is denoted as N_trace; the z-axis represents time, and the total number of time sampling points on a single trace is denoted as N_time.
[0058] Step S2: Set the current seismic trace line number as id_line and the seismic trace number as id_trace; extract the amplitude data of the current seismic trace from the three-dimensional seismic data volume and record it as a one-dimensional array x(n). The length of the one-dimensional array x(n) is equal to the total number of time sampling points N_time of the current seismic trace.
[0059] Step S3: Perform a -90° phase shift on the one-dimensional array x(n) to obtain the phase-shifted array y(n).
[0060] Specifically, the process of performing a -90° phase shift on the one-dimensional array x(n) is as follows:
[0061] Step S3.1: Perform a discrete Fourier transform on the one-dimensional array x(n) to obtain the spectrum array X(k), where K∈[0,N_time-1], and the part of k∈[(N_time-1) / 2,N_time-1] corresponds to the negative frequency.
[0062] Step S3.2: Perform frequency division calculations on the spectrum array X(k) to obtain the frequency-divided array Y(k), whose calculation formula is as follows:
[0063]
[0064] Step S3.3: Perform an inverse discrete Fourier transform on the frequency-division processed array Y(k). After the transform is completed, take its real part to obtain the phase-shifted array y(n).
[0065] Step S4: Perform time-frequency continuous wavelet transform on the phase shift array y(n) to obtain multiple common frequency amplitude arrays YAi(n) in the mid-to-high frequency band.
[0066] It should be noted that the frequency range of the above common frequency amplitude array YAi(n) is from 0.8 times the main frequency to 2 times the main frequency; among them, YA1(n) represents the amplitude array with a frequency of 0.8×Fre_main, YA2(n) represents the amplitude array with a frequency of 0.8×Fre_main+Δf, YA3(n) represents the amplitude array with a frequency of 0.8×Fre_main+2×Δf, ..., Δf takes the value 2 to 5.
[0067] Step S5: Perform Hilbert transform on each common frequency amplitude array YAi(n) to obtain the instantaneous phase array YA_Phasei(n) corresponding to each common frequency amplitude array YAi(n).
[0068] Specifically, the process of performing a Hilbert transform on the common frequency amplitude array YAi(n) is as follows:
[0069] Step S5.1: Perform a discrete Fourier transform on the common frequency amplitude array YAi(n) to obtain the spectrum array YA_Fi(k), where K∈[0,N_time-1] and the part of k∈[(N_time-1) / 2,N_time-1] corresponds to the negative frequency.
[0070] Step S5.2: Perform frequency division calculations on the spectrum array YA_Fi(k) to obtain the frequency-divided array YT. i (k).
[0071] Specifically, array YT i The formula for calculating (k) is:
[0072]
[0073] Step S5.3: Divide the frequency-division array YT i (k) Perform inverse discrete Fourier transform to obtain the analytic signal array YA_HTi(n). Record the real part data in the analytic signal array YA_HTi(n) as the real part array R_YA_HTi(n), and record the imaginary part data in the analytic signal array YA_HTi(n) as the imaginary part array I_YA_HTi(n).
[0074] Step S5.4: Calculate the instantaneous phase array YA_Phasei(n) based on the real part array R_YA_HTi(n) and the imaginary part array I_YA_HTi(n).
[0075] Specifically, the formula for calculating the instantaneous phase array YA_Phasei(n) is as follows:
[0076]
[0077] Step S6: For each instantaneous phase array YA_Phasei(n), calculate the second derivative Gi(n) with respect to the time variable n to obtain the three-dimensional data volume of the weak signal indication of the current seismic trace co-frequency.
[0078] Step S7: Change the seismic trace and repeat steps S2-S6 until all seismic traces in the study area are calculated, and obtain multiple weak signal indicators of the same frequency in the mid-to-high frequency band for three-dimensional data volume.
[0079] Step S8: Match all three-dimensional data volumes of weak signals with the same frequency with the top interface depth data and time-depth relationship data input in step S1. Select the three-dimensional data volume of weak signals with the same frequency that has the highest time-depth relationship with the weak reflection interface of the target section of the completed well in the study area. Record it as the optimal weak signal indicator data volume of the study area and complete the detection.
[0080] This invention is based on seismic data reconstruction using instantaneous phase. It establishes a correspondence between weak signals and characteristic data volumes by using techniques such as phase shifting, frequency division, data optimization, and second derivative calculation, based on the sensitivity of instantaneous phase to stratigraphic structure. This can effectively find three-dimensional data volumes that indicate weak signals and improve the accuracy of seismic weak signal detection.
[0081] Example 2
[0082] like Figure 2 As shown, this invention provides a seismic weak signal detection system based on instantaneous phase analysis, comprising:
[0083] The data input module is used to input the three-dimensional seismic data volume of the study area, to input the top interface depth data and time-depth relationship data of the target layer after drilling, and to record the dominant frequency of the seismic signal of the target layer in the three-dimensional seismic data volume as Fre_main.
[0084] The three-dimensional seismic data volume includes the x-axis, y-axis, and z-axis. The x-axis represents the direction of the seismic line, and the number of lines is denoted as N_line; the y-axis represents the direction of the seismic trace, and the total number of traces for a single seismic line is denoted as N_trace; the z-axis represents time, and the total number of time sampling points for a single seismic trace is denoted as N_time.
[0085] The setting and extraction module is used to set the current seismic trace line number as id_line and the seismic trace number as id_trace; it is used to extract the amplitude data of the current seismic trace and record it as a one-dimensional array x(n), the length of which is equal to the total number of time sampling points N_time of the current seismic trace.
[0086] The phase shift processing module is used to perform a -90° phase shift processing on the one-dimensional array x(n) and obtain the phase shift array y(n).
[0087] The transformation module is used to perform time-frequency continuous wavelet transform on the phase shift array y(n) and obtain multiple common frequency amplitude arrays YAi(n) in the mid-to-high frequency band.
[0088] The Hilbert transform module is used to perform Hilbert transform on each common frequency amplitude array YAi(n) and obtain the instantaneous phase array YA_Phasei(n) corresponding to each common frequency amplitude array YAi(n).
[0089] The second derivative calculation module is used to calculate the second derivative Gi(n) with respect to the time variable n for each instantaneous phase array YA_Phasei(n), and obtain the three-dimensional data volume of weak signal indication of the current seismic trace cofrequency.
[0090] The integrated processing module is used to replace seismic traces. It controls the setting and extraction module, phase shift processing module, transformation module, Hilbert transformation module, and second derivative calculation module to repeatedly process the replaced seismic traces until all seismic traces in the study area are calculated and a three-dimensional data volume of weak signals with multiple common frequencies in the mid-to-high frequency band is obtained.
[0091] The optimal matching module is used to match all three-dimensional data volumes of weak signals with the same frequency with the top interface depth data and time-depth relationship data in the data input module. It is used to select the three-dimensional data volume of weak signals with the same frequency as the weak reflection interface time-depth relationship of the target section of the completed well in the study area with the highest consistency rate, thereby obtaining the optimal weak signal indicator data volume of the study area.
[0092] In detail, each module in the system described in this embodiment uses the same technical means as in Embodiment 1 and can produce the same technical effects, which will not be repeated here.
[0093] Example 3
[0094] This invention provides an electronic device comprising:
[0095] At least one processor; and,
[0096] A memory communicatively connected to the at least one processor; wherein,
[0097] The memory stores a computer program that can be executed by the at least one processor, which enables the at least one processor to perform the aforementioned method for detecting weak seismic signals based on instantaneous phase analysis.
[0098] The present invention also provides a storage medium storing instructions that, when executed on a computer, cause the computer to perform the aforementioned method for detecting weak seismic signals based on instantaneous phase analysis.
[0099] Example 4
[0100] This embodiment uses a research area in Sichuan as an example to verify the detection method of Example 1, as detailed below:
[0101] Figure 3 The original seismic profiles of two completed wells (PY1 and PY3) in the study area are shown. The gray layers in the figures correctly identify the seismic horizons of the top interfaces of the Maokou Formation and the Maokou-2 Member. The seismic reflections from the top interface of the Maokou Formation exhibit strong peak reflection characteristics, making them easy to identify accurately, while the seismic reflections from the top interface of the Maokou-2 Member show typical weak reflection characteristics. Figure 3 The seismic horizon of the Maoerding interface could not be accurately traced; the short horizontal lines next to the two completed wells marked the longitudinal position of the actually drilled Maoerding interface.
[0102] Figure 4 It shows the corresponding Figure 3 The instantaneous phase profile, calculated using conventional methods, shows gray layers indicating the correct seismic horizons at the Maokouding and Maoerding interfaces. Short horizontal lines next to the two completed wells indicate the longitudinal position of the actually drilled Maoerding interface. However, from... Figure 4 It is also impossible to accurately track the seismic horizon at the top of Mao Erding.
[0103] Figure 5 The corresponding result is obtained using the detection method of this invention. Figure 3 The optimal weak signal indication profile at a frequency of 40Hz is shown in the image. The gray layers in the image correctly identify the seismic horizons of the top interface of the Maokou Formation and the top interface of the Maokou-2 Member. The short horizontal lines next to the two completed wells indicate the longitudinal position of the drilled Maokou-2 Member top interface. Figure 5 It can accurately track the seismic horizon at the top of Mao Erding.
[0104] As can be seen from this embodiment, by establishing a correspondence between weak signals and characteristic data volumes, the present invention can effectively identify and indicate signal reflections. Accordingly, the present invention provides technical support and basis for seismic horizon interpretation or reservoir prediction.
[0105] The above description is merely a specific embodiment of the present invention. Any feature disclosed in this specification may be replaced by other equivalent or similar features unless otherwise specified. All features or steps in the disclosed methods or processes may be combined in any way, except for mutually exclusive features and / or steps.
Claims
1. A seismic weak signal detection method based on instantaneous phase analysis, comprising the following steps: Step S1: Input the 3D seismic data volume of the study area, as well as the top interface depth data and time-depth relationship data of the target layer after drilling, and denote the dominant frequency of the seismic signal of the target layer in the 3D seismic data volume as Fre_main; wherein, the 3D seismic data volume includes the x-axis, y-axis and z-axis, the x-axis represents the direction of the survey line, and the number of lines of the survey line is denoted as N_line; the y-axis represents the direction of the seismic trace, and the total number of seismic traces of a single survey line is denoted as N_trace; the z-axis represents time, and the total number of time sampling points of a single seismic trace is denoted as N_time; Step S2: Set the current seismic trace line number as id_line and the seismic trace number as id_trace; extract the amplitude data of the current seismic trace from the three-dimensional seismic data volume and record it as a one-dimensional array x(n), the length of which is equal to the total number of time sampling points N_time of the current seismic trace; Step S3: Perform a -90° phase shift on the one-dimensional array x(n) to obtain the phase-shifted array y(n); Step S4: Perform time-frequency continuous wavelet transform on the phase shift array y(n) to obtain multiple common frequency amplitude arrays YAi(n) in the mid-to-high frequency band; Step S5: Perform Hilbert transform on each common frequency amplitude array YAi(n) to obtain the instantaneous phase array YA_Phasei(n) corresponding to each common frequency amplitude array YAi(n); Step S6: For each instantaneous phase array YA_Phasei(n), calculate the second derivative Gi(n) with respect to the time variable n to obtain the three-dimensional data volume of the weak signal indication of the current seismic trace co-frequency. Step S7: Change the seismic trace and repeat steps S2-S6 until all seismic traces in the study area are calculated, and obtain multiple weak signal indicators of the same frequency in the mid-to-high frequency band for three-dimensional data volume. Step S8: Match all three-dimensional data volumes of weak signals with the same frequency with the top interface depth data and time-depth relationship data input in step S1. Select the three-dimensional data volume of weak signals with the same frequency that has the highest time-depth relationship with the weak reflection interface of the target section of the completed well in the study area. Record it as the optimal weak signal indicator data volume of the study area and complete the detection.
2. The method for detecting weak seismic signals based on instantaneous phase analysis according to claim 1, characterized in that: In step S3, the process of performing a -90° phase shift on the one-dimensional array x(n) is as follows: Step S3.1: Perform a discrete Fourier transform on the one-dimensional array x(n) to obtain the spectrum array X(k), where K∈[0,N_time-1], and the part of k∈[(N_time-1) / 2,N_time-1] corresponds to the negative frequency; Step S3.2: Perform frequency division calculation on the spectrum array X(k) to obtain the frequency-divided array Y(k); Step S3.3: Perform an inverse discrete Fourier transform on the frequency-division processed array Y(k), take its real part, and obtain the phase-shifted array y(n).
3. The method for detecting weak seismic signals based on instantaneous phase analysis according to claim 2, characterized in that: In step S3.2, the formula for calculating the frequency division of the spectrum array X(k) is as follows:
4. The method for detecting weak seismic signals based on instantaneous phase analysis according to claim 3, characterized in that: In step S4, the frequency range of the common frequency amplitude array YAi(n) is from 0.8 times the main frequency to 2 times the main frequency; where YA1(n) represents the amplitude array with a frequency of 0.8×Fre_main, YA2(n) represents the amplitude array with a frequency of 0.8×Fre_main+Δf, YA3(n) represents the amplitude array with a frequency of 0.8×Fre_main+2×Δf, ..., Δf takes the value 2 to 5.
5. The method for detecting weak seismic signals based on instantaneous phase analysis according to claim 4, characterized in that: In step S5, the process of performing a Hilbert transform on the common frequency amplitude array YAi(n) is as follows: Step S5.1: Perform a discrete Fourier transform on the common frequency amplitude array YAi(n) to obtain the spectrum array YA_Fi(k), where K∈[0,N_time-1], and the part of k∈[(N_time-1) / 2,N_time-1] corresponds to the negative frequency; Step S5.2, frequency division calculation is performed on the spectrum array YA_Fi(k) to obtain a frequency-division processed array YT i (k); Step S5.3, processing the array YT after frequency division i (k) performing inverse discrete Fourier transform to obtain an analytic signal array YA_HTi(n), taking real part data in the analytic signal array YA_HTi(n) as a real part array R_YA_HTi(n), and taking imaginary part data in the analytic signal array YA_HTi(n) as an imaginary part array I_YA_HTi(n); Step S5.4: Calculate the instantaneous phase array YA_Phasei(n) based on the real part array R_YA_HTi(n) and the imaginary part array I_YA_HTi(n).
6. The method for detecting weak seismic signals based on instantaneous phase analysis according to claim 5, characterized in that: In step S5.2, the formula for calculating the frequency division of the spectrum array YA_Fi(k) is as follows:
7. The method for detecting weak seismic signals based on instantaneous phase analysis according to claim 5, characterized in that: In step S5.4, the formula for calculating the instantaneous phase array YA_Phasei(n) is as follows:
8. A system for detecting seismic weak signals based on instantaneous phase analysis, characterized in that include: The data input module is used to input the three-dimensional seismic data volume of the study area, to input the top interface depth data and time-depth relationship data of the target layer after drilling, and to record the dominant frequency of the seismic signal of the target layer in the three-dimensional seismic data volume as Fre_main; The setting and extraction module is used to set the current seismic trace's survey line number as id_line and the seismic trace number as id_trace; it is used to extract the amplitude data of the current seismic trace and record it as a one-dimensional array x(n), the length of which is equal to the total number of time sampling points N_time of the current seismic trace; The phase shift processing module is used to perform a -90° phase shift processing on the one-dimensional array x(n) and obtain the phase shift array y(n); The transformation module is used to perform time-frequency continuous wavelet transform on the phase shift array y(n) and obtain multiple common frequency amplitude arrays YAi(n) in the mid-to-high frequency band. The Hilbert transform module is used to perform Hilbert transform on each common frequency amplitude array YAi(n) and obtain the instantaneous phase array YA_Phasei(n) corresponding to each common frequency amplitude array YAi(n); The second derivative calculation module is used to calculate the second derivative Gi(n) with respect to the time variable n for each instantaneous phase array YA_Phasei(n), and obtain the three-dimensional data volume of weak signal indication of the current seismic trace co-frequency; The integrated processing module is used to replace the seismic traces. It controls the setting and extraction module, phase shift processing module, transformation module, Hilbert transformation module and second derivative calculation module to repeatedly process the replaced seismic traces until all seismic traces in the study area are calculated and multiple weak signal indicators in the mid-to-high frequency band are obtained. The optimal matching module is used to match all three-dimensional data volumes of weak signals with the same frequency with the top interface depth data and time-depth relationship data in the data input module. It is used to select the three-dimensional data volume of weak signals with the same frequency as the weak reflection interface time-depth relationship of the target section of the completed well in the study area with the highest consistency rate, thereby obtaining the optimal weak signal indicator data volume of the study area.
9. An electronic device, comprising: The electronic device includes: At least one processor; and, A memory communicatively connected to the at least one processor; wherein, The memory stores a computer program that can be executed by the at least one processor, the computer program being executed by the at least one processor to enable the at least one processor to perform the seismic weak signal detection method based on instantaneous phase analysis as described in any one of claims 1-7.
10. A storage medium characterized by: The storage medium stores instructions that, when executed on a computer, cause the computer to perform the seismic weak signal detection method based on instantaneous phase analysis as described in any one of claims 1-7.