A method for enhancing weak seismic signals

Abstract

The present application relates to the technical field of seismic data processing, and specifically discloses a seismic weak signal enhancement processing method, which utilizes the physical characteristics that effective weak signals are continuous and consistent in space phase, while random noise phase randomly jumps, performs coherent stacking along a track, so that the effective weak signals obtain positive energy compensation, while the noise cannot obtain gain, thereby passively suppressing the noise in a relative sense, significantly improving the signal-to-noise ratio and fidelity of deep weak signals, constructing a time-frequency topology network based on a phase locking value by extracting an instantaneous phase, and extracting a spatial track in combination with a graph search algorithm, thereby breaking through the bottleneck of poor adaptability of traditional methods to complex geological structures.

Description

Technical Field

[0001] This invention relates to the field of earthquake data processing technology, and in particular to a method for enhancing weak earthquake signals. Background Technology

[0002] Currently, seismic waves undergo strong absorption and attenuation in deep oil and gas exploration, resulting in extremely weak effective reflected signal energy returning to the surface. These signals are often submerged by background random noise or shielded by strong shallow reflections. Accurately extracting and enhancing these weak deep signals is a prerequisite for improving the imaging accuracy of complex underground structures and the accuracy of reservoir characterization.

[0003] However, existing technologies mainly employ passive separation methods such as transform domain threshold filtering, strong reflection wavelet stripping, or matrix rank reduction, treating these methods as separation objects based solely on energy differences or horizontal linear characteristics. This completely ignores the physical continuity and phase consistency of weak signals extending along space in complex structures. Furthermore, since strong layer stripping methods use subtraction operations, they can erase the real weak signals immediately below due to wavelet superposition effects. Moreover, matrix rank reduction is highly dependent on the assumption of a horizontal phase axis, making it impossible to actively track the spatial trajectory of the weak signal's phase axis and achieve directional energy compensation in complex structural regions without damaging the dynamic characteristics of the weak signal.

[0004] Therefore, there is an urgent need for a method to enhance weak seismic signals in order to solve the problems in the existing technologies mentioned above. Summary of the Invention

[0005] The purpose of this invention is to provide a method for enhancing weak seismic signals, aiming to solve the technical problems in the prior art where weak signals are easily mistakenly ignored, strong layer stripping easily causes collateral damage, and the adaptability to complex structures is poor.

[0006] To achieve the above objectives, the present invention employs a seismic weak signal enhancement processing method, comprising the following steps: Time-frequency transformation is performed on noisy seismic data to obtain the instantaneous amplitude matrix and instantaneous phase matrix of the time spectrum; Calculate the instantaneous phase difference between adjacent seismic traces at the same time and frequency nodes, and construct a time-frequency topology network characterizing spatial phase consistency based on the instantaneous phase difference; In the time-frequency topology network, strong reflection nodes are eliminated based on the local root mean square amplitude adaptive threshold. In the remaining low-amplitude subnetwork, a graph search algorithm is used to extract the spatial trajectory of the weak reflection in-phase axis with continuous phase consistency. Along the spatial trajectory of the weak reflection in-phase axis, the nodes in the neighborhood are weighted coherently superimposed based on the phase coherence factor to obtain the enhanced time spectrum; The enhanced time spectrum is subjected to inverse time-frequency transformation to reconstruct the weak signal enhanced seismic data.

[0007] The time-frequency transformation employs a generalized S-transform, and an adaptive time-frequency window that varies non-linearly with frequency is introduced during the transformation process. The width of the adaptive time-frequency window decreases as the frequency increases, in order to balance the time resolution of low-frequency weak signals with the frequency resolution of high-frequency weak signals.

[0008] The calculation of the instantaneous phase difference and the quantization of spatial phase consistency employ a phase-locked value, the calculation formula of which is as follows: In the formula, Adjacent seismic traces and The phase-locked value at time t and frequency f, where W is the length of the time sliding window. The time sampling point number within the sliding window. For instantaneous phase, The imaginary unit; When adjacent channels are in phase, the phase lock value approaches 1; when random noise or faults are present, the phase lock value approaches 0.

[0009] The time-frequency topology network is an undirected weighted graph structure, where the vertices of the undirected weighted graph are data points in the time-frequency spectrum, and the weights of the edges of the undirected weighted graph are determined by the product of the phase-locking value between adjacent channels and the instantaneous amplitude similarity.

[0010] The method for calculating the adaptive threshold of the local root mean square amplitude is as follows: calculate the local moving root mean square value of the amplitude of the time-frequency node along the spatial direction, and determine the node whose instantaneous amplitude is greater than 1.5 to 2.0 times the local moving root mean square value as a strong reflection node and remove it from the time-frequency topology network.

[0011] The specific process of extracting the spatial trajectory of the weak reflection phase axis using the graph search algorithm is as follows: set the phase lock value threshold to 0.7~0.8, retain the edges with the phase lock value greater than the threshold in the low amplitude sub-network; use the depth-first search algorithm to traverse the connected components of the remaining network, filter out short branches with a spatial extension length less than 3 consecutive seismic traces, and define the retained long-range connected components as the spatial trajectory of the weak reflection phase axis.

[0012] The calculation formula for the weighted coherent superposition based on the phase coherence factor is as follows: In the formula, For the enhanced time-frequency amplitude, The original time-frequency amplitude, For adaptive gain coefficients, Let be the set of neighborhoods of the weakly reflecting in-phase axis spatial trajectory of the current node. For the node index within the neighborhood, For Gaussian decreasing weights based on spatial distance, This represents the instantaneous phase difference between neighboring nodes and the current node. This is the phase coherence factor.

[0013] The phase coherence factor is used to enable effective weak reflection signals with consistent phase to obtain positive superposition gain. When random noise exists in the neighborhood, the instantaneous phase difference of the random noise is randomly distributed, causing the positive and negative values ​​of the phase coherence factor to cancel each other out, resulting in the random noise amplitude being unable to obtain superposition gain, thereby achieving passive suppression of random noise.

[0014] In the process of performing the weighted coherent superposition, a nonlinear saturation limiting function is introduced to limit the maximum amplitude gain of a single superposition. The limiting threshold of the nonlinear saturation limiting function is dynamically determined based on the difference between the local noise variance and the local signal variance, so as to prevent the waveform distortion or the generation of false in-phase axes due to overcompensation of weak signal amplitude.

[0015] Among them, time-frequency nodes that are not covered by any weak reflection phase axis spatial trajectory retain their original amplitude and do not participate in weighted coherent superposition; after processing all time-frequency nodes, the processed time spectrum is reconstructed into a spatiotemporal seismic profile using generalized S-inverse transform.

[0016] The earthquake weak signal enhancement processing method of the present invention has the following beneficial effects: 1. This invention utilizes the physical characteristics that effective weak signals have continuous and consistent phase in space, while random noise has randomly jumping phase. Coherent superposition is performed along the trajectory, so that the effective weak signals can obtain positive energy compensation, while the noise cannot obtain gain. Thus, the noise is passively suppressed in a relative sense, and the signal-to-noise ratio and fidelity of deep weak signals are significantly improved. 2. Instead of fitting and subtracting the strong reflection layer, this invention directly removes high-energy nodes from the topology network and actively searches for weak signal trajectories in the remaining low-amplitude pure noise background. This cuts off the physical damage path of strong layer processing to weak layer signals and protects the real weak reflection waveform to the greatest extent. 3. This invention constructs a time-frequency topology network based on phase-locked values ​​by extracting instantaneous phases and extracts spatial trajectories by combining graph search algorithms. This overcomes the bottleneck of poor adaptability of traditional methods to complex geological structures. The phase-locked values ​​constructed in this invention measure the dynamic continuity of wave field propagation between adjacent channels, rather than geometric horizontality. Even if the in-phase axis undergoes drastic tilt changes or breaks, as long as the wave field is physically continuous, a connected component can be formed in the topology network, thereby achieving accurate tracking and enhancement of weak signal in-phase axes under any complex structure. Detailed Implementation

[0017] This invention provides a method for enhancing weak seismic signals, comprising the following steps: Time-frequency transformation is performed on noisy seismic data to obtain the instantaneous amplitude matrix and instantaneous phase matrix of the time spectrum; Calculate the instantaneous phase difference between adjacent seismic traces at the same time and frequency nodes, and construct a time-frequency topology network representing spatial phase consistency based on the instantaneous phase difference; In the time-frequency topology network, strong reflection nodes are eliminated based on the local root mean square amplitude adaptive threshold. In the remaining low-amplitude subnetwork, a graph search algorithm is used to extract the spatial trajectory of the weak reflection in-phase axis with continuous phase consistency. Along the spatial trajectory of the weak reflection in-phase axis, the nodes in the neighborhood are weighted coherently superimposed based on the phase coherence factor to obtain the enhanced time spectrum; The enhanced time spectrum is subjected to inverse time-frequency transformation to reconstruct the weak signal enhanced seismic data.

[0018] Furthermore, the time-frequency transformation adopts the generalized S-transform, and an adaptive time-frequency window that changes nonlinearly with frequency is introduced during the transformation process. The width of the adaptive time-frequency window decreases as the frequency increases, so as to balance the time resolution of low-frequency weak signals and the frequency resolution of high-frequency weak signals.

[0019] Furthermore, the calculation of the instantaneous phase difference and the quantization of spatial phase consistency adopt a phase-locked value, the calculation formula of which is: In the formula, Adjacent seismic traces and The phase-locked value at time t and frequency f, where W is the length of the time sliding window. The time sampling point number within the sliding window. For instantaneous phase, The imaginary unit; When adjacent channels are in phase, the phase lock value approaches 1; when random noise or faults are present, the phase lock value approaches 0.

[0020] Furthermore, the time-frequency topology network is an undirected weighted graph structure, where the vertices of the undirected weighted graph are data points in the time-frequency spectrum, and the weights of the edges of the undirected weighted graph are determined by the product of the phase-locking value between adjacent channels and the instantaneous amplitude similarity.

[0021] Furthermore, the method for calculating the local root mean square amplitude adaptive threshold is as follows: calculate the local moving root mean square value of the amplitude of the time-frequency node along the spatial direction, and determine the node whose instantaneous amplitude is greater than 1.5 to 2.0 times the local moving root mean square value as a strong reflection node and remove it from the time-frequency topology network.

[0022] Furthermore, the specific process of extracting the spatial trajectory of the weak reflection phase axis using the graph search algorithm is as follows: set the phase lock value threshold to 0.7~0.8, retain the edges with the phase lock value greater than the threshold in the low amplitude sub-network; use the depth-first search algorithm to traverse the connected components of the remaining network, filter out short branches with a spatial extension length less than 3 consecutive seismic traces, and define the retained long-range connected components as the spatial trajectory of the weak reflection phase axis.

[0023] Furthermore, the calculation formula for the weighted coherent superposition based on the phase coherence factor is as follows: In the formula, For the enhanced time-frequency amplitude, The original time-frequency amplitude, For adaptive gain coefficients, Let be the set of neighborhoods of the weakly reflecting in-phase axis spatial trajectory of the current node. For the node index within the neighborhood, For Gaussian decreasing weights based on spatial distance, This represents the instantaneous phase difference between neighboring nodes and the current node. This is the phase coherence factor.

[0024] Furthermore, the phase coherence factor is used to enable effective weak reflection signals with consistent phase to obtain positive superposition gain; when random noise exists in the neighborhood, the instantaneous phase difference of the random noise is randomly distributed, causing the positive and negative values ​​of the phase coherence factor to cancel each other out, resulting in the random noise amplitude being unable to obtain superposition gain, thereby achieving passive suppression of random noise.

[0025] Furthermore, when performing the weighted coherent superposition, a nonlinear saturation limiting function is introduced to limit the maximum amplitude gain of a single superposition. The limiting threshold of the nonlinear saturation limiting function is dynamically determined based on the difference between the local noise variance and the local signal variance, in order to prevent the weak signal amplitude from being overcompensated, which could lead to waveform distortion or the generation of false in-phase axes.

[0026] Furthermore, time-frequency nodes not covered by any weakly reflecting phase axis spatial trajectories retain their original amplitudes and do not participate in weighted coherent superposition; after processing all time-frequency nodes, the processed time-frequency spectrum is reconstructed into a spatiotemporal seismic profile using a generalized S-transform.

[0027] In this embodiment, This invention provides a method for enhancing weak seismic signals, comprising the following steps: S1: Perform time-frequency transformation on noisy seismic data to obtain the instantaneous amplitude matrix and instantaneous phase matrix of the time spectrum; Furthermore, the time-frequency transformation adopts the generalized S-transform, and an adaptive time-frequency window that changes nonlinearly with frequency is introduced during the transformation process. The width of the adaptive time-frequency window decreases as the frequency increases, so as to balance the time resolution of low-frequency weak signals and the frequency resolution of high-frequency weak signals. Given noisy seismic data Where s is the spatial channel number and t is time, the time spectrum is obtained using the generalized S-transform. Then, the instantaneous amplitude matrix is ​​obtained by analyzing the signal. and instantaneous phase matrix .

[0028] S2: Calculate the instantaneous phase difference between adjacent seismic traces at the same time and frequency nodes, and construct a time-frequency topology network representing spatial phase consistency based on the instantaneous phase difference; Furthermore, the calculation of the instantaneous phase difference and the quantization of spatial phase consistency adopt a phase-locked value, the calculation formula of which is: In the formula, Adjacent seismic traces and The phase-locked value at time t and frequency f, where W is the length of the time sliding window. The time sampling point number within the sliding window. For instantaneous phase, The value is an imaginary unit; when adjacent channels are in phase, the phase lock value approaches 1, and when random noise or faults are present, the phase lock value approaches 0. It should be noted that the time variable t and the sliding window length W involved in the formula are both discrete integers. When W is odd, the upper and lower limits of the summation are... The floor function is used to ensure that the summation interval contains valid integer indices. Furthermore, the time-frequency topology network is an undirected weighted graph structure. The vertex V of the undirected weighted graph is a data point in the time spectrum, and the weight of the edge E of the undirected weighted graph is determined by the product of the phase lock value between adjacent channels and the instantaneous amplitude similarity. The instantaneous amplitude similarity S is calculated using Gaussian radial basis functions: In the formula, For local amplitude variance, the final weight of the edge This simultaneously constrains the spatial phase continuity and waveform amplitude similarity of weak signals.

[0029] S3: In the time-frequency topology network, strong reflection nodes are eliminated based on the local root mean square amplitude adaptive threshold. In the remaining low amplitude subnetwork, the spatial trajectory of weak reflection in phase axis with continuous phase consistency is extracted using a graph search algorithm. Furthermore, the method for calculating the local root-mean-square amplitude adaptive threshold is as follows: calculate the local sliding root-mean-square value of the time-frequency node amplitude along the spatial direction. In the formula, M represents the number of channels in the spatial sliding window, which represents the instantaneous amplitude. Nodes with values ​​greater than 1.5 to 2.0 times the local moving root mean square value are identified as strong reflection nodes and removed from the time-frequency topology network to avoid subsequent subtraction stripping causing collateral damage to the underlying weak signal. Furthermore, the specific process of extracting the spatial trajectory of the weak reflection phase axis using the graph search algorithm is as follows: set the phase lock value threshold to 0.7~0.8, and retain the edges with the phase lock value greater than the threshold in the low amplitude sub-network; use the DFS algorithm to traverse the connected components of the remaining network. In the DFS traversal, the search direction is only along the spatial direction (i.e., adjacent seismic traces), and find connected nodes with the same frequency and time shift in the time direction less than half a phase period between adjacent traces as the next traversal nodes. Filter out short branches with a spatial extension length less than 3 consecutive seismic traces, and define the retained long-range connected components as the spatial trajectory of the weak reflection phase axis.

[0030] S4: Along the spatial trajectory of the weak reflection in-phase axis, perform weighted coherent superposition of the nodes in the neighborhood based on the phase coherence factor to obtain the enhanced time spectrum; Furthermore, the calculation formula for the weighted coherent superposition based on the phase coherence factor is as follows: In the formula, For the enhanced time-frequency amplitude, The original time-frequency amplitude, For adaptive gain coefficients, Let be the set of neighborhoods of the weakly reflecting in-phase axis spatial trajectory of the current node. For the node index within the neighborhood, For Gaussian decreasing weights based on spatial distance, This represents the instantaneous phase difference between neighboring nodes and the current node. The phase coherence factor; In specific implementation, the neighborhood set Defined as the set of nodes within N consecutive spatial channels (e.g., N=5) along the spatial trajectory of the same phase axis, the adaptive gain coefficient... Based on the dynamic determination of the local signal-to-noise ratio estimate, the formula is as follows: ,in, The basic gain constant, To prevent division by zero, a regularization constant is used; Furthermore, the phase coherence factor is used to enable effective weak reflection signals with consistent phase to obtain positive superposition gain; when random noise exists in the neighborhood, the instantaneous phase difference of the random noise is randomly distributed, causing the positive and negative values ​​of the phase coherence factor to cancel each other out, resulting in the random noise amplitude being unable to obtain superposition gain, thereby achieving passive suppression of random noise. Furthermore, during the weighted coherent superposition, a nonlinear saturation limiting function is introduced to restrict the maximum amplitude gain of a single superposition. The limiting threshold of the nonlinear saturation limiting function is dynamically determined based on the difference between the local noise variance and the local signal variance to prevent overcompensation of weak signal amplitude, which could lead to waveform distortion or the generation of false in-phase axes. Specifically, the local noise variance is estimated using the robust median absolute deviation method, with the following formula: Local signal variance through We obtain, where The total variance of the local data, and the amplitude limit threshold. When the compensated amplitude Exceed At that time, it was forcibly truncated into .

[0031] S5: Perform inverse time-frequency transformation on the enhanced time spectrum to reconstruct the weak signal enhanced seismic data; Furthermore, time-frequency nodes not covered by any weakly reflecting phase axis spatial trajectories retain their original amplitudes and do not participate in weighted coherent superposition; after processing all time-frequency nodes, the processed time-frequency spectrum is reconstructed into a spatiotemporal seismic profile using a generalized S-transform.

[0032] I. Time-Frequency Decomposition Module (Feature Extraction Layer) 1. Function and scenario integration: In seismic exploration, deep weak signals often manifest as non-stationary transient signals. This module is used to convert one-dimensional seismic trace data in the spatiotemporal domain into a two-dimensional complex matrix in the time and frequency domain, so as to extract instantaneous phase features that are not affected by energy attenuation and instantaneous amplitude features that characterize energy.

[0033] 2. Algorithm details: Noisy seismic data is processed using a generalized S-transform with an adaptive time-frequency window. For a given spatiotemporal domain seismic data... Where x is the spatial channel number and t is time, its generalized S-transform formula is: In the formula, f is the frequency. For time integration variables, An adaptive window function that varies non-linearly with frequency f has a window width that decreases as frequency increases in order to balance the time resolution of low-frequency weak signals with the frequency resolution of high-frequency weak signals. The instantaneous amplitude matrix of the time spectrum is extracted by analyzing the signal. and instantaneous phase matrix .

[0034] 3. Input / Output Input data: Noisy seismic data in the spatiotemporal domain .

[0035] Output data: Instantaneous amplitude matrix in the time-frequency domain and instantaneous phase matrix .

[0036] II. Topology Network Construction Module (Topology Mapping Layer) 1. Function and scenario integration: This module maps discrete data points in the time and frequency domain to an undirected weighted graph network in graph theory. In the earthquake physics scenario, reflected waves from the same stratum interface have a constant phase difference (i.e., continuous phase axis) when they reach adjacent detectors. This module constructs a topological network by quantifying this spatial phase consistency.

[0037] 2. Algorithm details and model structure: Construct an undirected weighted graph ; Vertex set V: Each data point in the time spectrum Defined as a vertex.

[0038] Edge set E and weight calculation: Only edges between adjacent seismic traces (x and x+1) at the same frequency f and similar time t are calculated, and the edge weights are calculated. It is determined by the product of the phase consistency factor and the amplitude similarity factor.

[0039] (1) Phase Lock Value (PLV) Algorithm: Used to quantify spatial phase consistency, the formula is: In the formula, W is the length of the time sliding window (discrete integer), n is the index of the time sampling point within the sliding window, and i is the imaginary unit. It should be noted that when W is odd, the upper and lower limits of the summation are... Using the floor operation, when adjacent channels are in phase and the phase difference is constant, the magnitude after summation approaches 1. When there is random noise or fault, the phase difference changes randomly, and the magnitude approaches 0.

[0040] (2) Instantaneous amplitude similarity algorithm: The Gaussian radial basis function is used for calculation, and the formula is: In the formula, This formula constrains the similarity of waveform amplitudes by representing the local amplitude variance.

[0041] (3) Weight fusion: .

[0042] 3. Input / Output: Input data: Instantaneous amplitude matrix and instantaneous phase matrix .

[0043] Output data: Time-frequency topological network (undirected weighted graph) characterizing the physical continuity of the space wave field. ).

[0044] III. Trajectory Tracking Module (Logic Processing Layer) 1. Function and scenario integration: This module aims to isolate strong reflection shielding layers and actively track the physical spatial trajectory of weak signal in-phase axes in pure noise background. Traditional methods rely on horizontal assumptions, while this algorithm adapts to steep structures through graph search.

[0045] 2. Algorithm details and model structure It contains two cascaded submodules: (1) Strong Layer Bypass Algorithm Submodule: Calculates the local moving root mean square (RMS) value of the time-frequency node amplitude along the spatial direction: In the formula, M is the number of channels in the spatial sliding window, which represents the instantaneous amplitude. Nodes with amplitudes greater than 1.5 to 2.0 times the RMS value are identified as strong reflection nodes and removed from graph G, generating a low-amplitude subgraph. This mechanism avoids the simultaneous erasure of underlying weak signals caused by traditional subtraction stripping.

[0046] (2) Depth-first search graph traversal algorithm submodule: in low amplitude subgraphs In this study, a PLV threshold of 0.7-0.8 is set, and only edges with weights greater than this threshold are retained. A Depth-First Search (DFS) algorithm is used to traverse connected components. The search rules for the seismic scenario are defined as follows: the search direction advances only along the spatial direction (trace number x increases); and between adjacent traces, the search seeks traces where the time shift in the temporal direction is less than half a phase period (i.e.,...). ,in Connected nodes with the same frequency (corresponding to the period of frequency f) are used as the nodes for the next traversal, in order to adapt to the time delay of the same phase axis caused by the tilt of the strata and filter out short branches (isolated noise) with a spatial extension length less than 3 consecutive seismic traces.

[0047] 3. Input / Output Input data: Time-frequency topology network Instantaneous amplitude matrix .

[0048] Output data: A set of spatial trajectories of weakly reflecting in-phase axes, i.e., an ordered sequence of three-dimensional coordinates. .

[0049] IV. Coherence Compensation and Reconstruction Module (Energy Compensation and Output Layer) 1. Function and scenario integration: This module performs directional energy superposition (addition) along the tracked physical trajectory, so that the effective weak signal can obtain positive gain, while random noise cannot obtain gain due to the positive and negative cancellation of the phase coherence factor, and is reconstructed back to the spatiotemporal domain after compensation.

[0050] 2. Algorithm details and model structure: It contains three sub-modules: (1) Weighted coherent superposition algorithm based on phase coherence factor: along the extracted trajectory, the nodes in the neighborhood are weighted and accumulated, and the formula is: In the formula, This is the set of indices of the neighboring nodes within N consecutive spatial channels before and after the current node along the trajectory. Gaussian decreasing weights based on spatial distance; Phase coherence factor, effective weak signal →1 Gain gain, random noise The positive and negative cancellation results in a gain of 0.

[0051] (2) Adaptive gain and limiting algorithm: To prevent waveform distortion caused by overcompensation of deep, extremely weak signals, an adaptive gain coefficient is introduced. Nonlinear saturation limiting, local noise variance Robust Median Absolute Deviation (MAD) estimation is used: Local signal variance through We obtain, where The total variance of the local data, and the amplitude limit threshold. When the compensated amplitude Exceed At that time, it was forcibly truncated into .

[0052] (3) Inverse time-frequency transform algorithm: For pure noise time-frequency nodes that are not covered by any trajectory, their original amplitude is kept and not superimposed. After the full-band nodes are processed, the generalized S-transform is used to reconstruct the time-frequency spectrum into spatiotemporal seismic data.

[0053] 3. Input / Output Input data: a set of spatial trajectories of weakly reflected in-phase axes, the original instantaneous amplitude matrix A, and the instantaneous phase matrix. .

[0054] Output data: Final seismic data volume after spatiotemporal domain weak signal enhancement .

[0055] The above description discloses only one preferred embodiment of the present invention, and should not be construed as limiting the scope of the present invention. Those skilled in the art will understand that all or part of the processes of the above embodiments can be implemented, and equivalent changes made in accordance with the claims of the present invention are still within the scope of the invention.

Claims

1. A method for enhancing weak seismic signals, characterized in that, Includes the following steps: Time-frequency transformation is performed on noisy seismic data to obtain the instantaneous amplitude matrix and instantaneous phase matrix of the time spectrum; Calculate the instantaneous phase difference between adjacent seismic traces at the same time and frequency nodes, and construct a time-frequency topology network representing spatial phase consistency based on the instantaneous phase difference; In the time-frequency topology network, strong reflection nodes are eliminated based on the local root mean square amplitude adaptive threshold. In the remaining low-amplitude subnetwork, a graph search algorithm is used to extract the spatial trajectory of the weak reflection in-phase axis with continuous phase consistency. Along the spatial trajectory of the weak reflection in-phase axis, the nodes in the neighborhood are weighted coherently superimposed based on the phase coherence factor to obtain the enhanced time spectrum; The enhanced time spectrum is subjected to inverse time-frequency transformation to reconstruct the weak signal enhanced seismic data.

2. The seismic weak signal enhancement processing method as described in claim 1, characterized in that, The time-frequency transformation adopts the generalized S-transform, and an adaptive time-frequency window that changes nonlinearly with frequency is introduced during the transformation process. The width of the adaptive time-frequency window decreases as the frequency increases, so as to balance the time resolution of low-frequency weak signals and the frequency resolution of high-frequency weak signals.

3. The seismic weak signal enhancement processing method as described in claim 1, characterized in that, The calculation of the instantaneous phase difference and the quantization of spatial phase consistency adopt the phase-locked value, the calculation formula of which is: In the formula, Adjacent seismic traces and The phase-locked value at time t and frequency f, where W is the length of the time sliding window. The time sampling point number within the sliding window. For instantaneous phase, The imaginary unit; When adjacent channels are in phase, the phase lock value approaches 1; when random noise or faults are present, the phase lock value approaches 0.

4. The seismic weak signal enhancement processing method as described in claim 3, characterized in that, The time-frequency topology network is an undirected weighted graph structure. The vertices of the undirected weighted graph are data points in the time spectrum. The weights of the edges of the undirected weighted graph are determined by the product of the phase-locking value between adjacent channels and the instantaneous amplitude similarity.

5. The seismic weak signal enhancement processing method as described in claim 1, characterized in that, The method for calculating the adaptive threshold of the local root mean square amplitude is as follows: calculate the local moving root mean square value of the amplitude of the time-frequency node along the spatial direction, and determine the node whose instantaneous amplitude is greater than 1.5 to 2.0 times the local moving root mean square value as a strong reflection node and remove it from the time-frequency topology network.

6. The seismic weak signal enhancement processing method as described in claim 4, characterized in that, The specific process of extracting the spatial trajectory of the weak reflection phase axis using the graph search algorithm is as follows: set the phase lock value threshold to 0.7~0.8, and retain the edges with the phase lock value greater than the threshold in the low amplitude sub-network; use the depth-first search algorithm to traverse the connected components of the remaining network, filter out the short branches with a spatial extension length less than 3 consecutive seismic traces, and define the retained long-range connected components as the spatial trajectory of the weak reflection phase axis.

7. The seismic weak signal enhancement processing method as described in claim 6, characterized in that, The calculation formula for the weighted coherent superposition based on the phase coherence factor is as follows: In the formula, For the enhanced time-frequency amplitude, The original time-frequency amplitude, For adaptive gain coefficients, Let be the set of neighborhoods of the weakly reflecting in-phase axis spatial trajectory of the current node. For the node index within the neighborhood, For Gaussian decreasing weights based on spatial distance, This represents the instantaneous phase difference between neighboring nodes and the current node. This is the phase coherence factor.

8. The seismic weak signal enhancement processing method as described in claim 7, characterized in that, The phase coherence factor is used to enable effective weak reflection signals with consistent phase to obtain positive superposition gain. When random noise exists in the neighborhood, the instantaneous phase difference of the random noise is randomly distributed, causing the positive and negative values ​​of the phase coherence factor to cancel each other out, resulting in the random noise amplitude not being able to obtain superposition gain, thereby achieving passive suppression of random noise.

9. The seismic weak signal enhancement processing method as described in claim 7, characterized in that, When performing the weighted coherent superposition, a nonlinear saturation limiting function is introduced to limit the maximum amplitude gain of a single superposition. The limiting threshold of the nonlinear saturation limiting function is dynamically determined based on the difference between the local noise variance and the local signal variance to prevent overcompensation of weak signal amplitude from causing waveform distortion or generating false in-phase axes.

10. The seismic weak signal enhancement processing method as described in claim 1, characterized in that, Time-frequency nodes not covered by any weakly reflecting phase axis spatial trajectories retain their original amplitudes and do not participate in weighted coherent superposition. After processing all time-frequency nodes, the processed time-frequency spectrum is reconstructed into a spatiotemporal seismic profile using a generalized S-transform.