Seismic interference signal enhancement method and system based on adaptive stacking strategy

Abstract

The application belongs to the technical field of seismic interference signal processing, and is a seismic interference signal enhancement method and system based on an adaptive stacking strategy, which comprises the following steps: converting a seismic interference time-domain signal to a frequency domain, and calculating a cross-correlation function data segment between each station pair; screening the cross-correlation function data segment; transforming the screened multi-channel cross-correlation function data segment to a frequency domain, and constructing a cross-spectrum matrix at each frequency point; designing an adaptive filter according to the cross-spectrum matrix; filtering the cross-correlation function data segment by using the adaptive filter to obtain a filtered signal; performing time-frequency decomposition on the filtered signal, and extracting a time transient phase; constructing a time-frequency domain coherence weight based on the transient phase information; and weighting and stacking the filtered signal by using the time-frequency domain coherence weight to obtain a stacked time-frequency spectrum. The application can effectively reduce noise interference in the signal, and make the signal waveform more stable and regular.

Description

Technical Field

[0001] This application belongs to the field of seismic interferometric signal processing technology, specifically a seismic interferometric signal enhancement method and system based on an adaptive superposition strategy. Background Technology

[0002] Seismic interferometry, which reconstructs empirical Green's functions from background noise cross-correlation functions, has become a key tool for detecting the fine structure of the Earth's crust and monitoring seismic activity. Traditional methods rely on the superposition of continuous noise data over months or even years to obtain cross-correlation functions with high signal-to-noise ratios (SNRs) and extract reliable surface wave dispersion information. However, in actual observation environments, background noise sources are often non-uniformly distributed, and factors such as station layout, observation period, and environmental interference limit the accumulation of coherent signals. Steady-state scattered signals in the cross-correlation function are often obscured by non-steady-state noise, significantly reducing the SNR. Although existing studies have employed linear superposition or time-frequency domain phase-weighted superposition to improve signal quality, traditional linear superposition has limited ability to suppress incoherent noise. Furthermore, conventional phase-weighted superposition uses fixed parameter settings, making it difficult to distinguish between effective wavefields and interference when noise sources are unevenly distributed or observation time is insufficient. This leads to jumps or decreased resolution in the final extracted dispersion curve, limiting the application of this technique in short-term pre-earthquake monitoring or rapid structural imaging. Summary of the Invention

[0003] To address the problems of low signal-to-noise ratio and poor stability of cross-correlation functions caused by uneven noise source distribution and insufficient observation time in existing technologies, this application provides a seismic interferometric signal enhancement method and system based on an adaptive superposition strategy.

[0004] A seismic interferometric signal enhancement method based on an adaptive superposition strategy, according to an embodiment of the first aspect of this application, includes:

[0005] The seismic interferometric time-domain signal was converted to the frequency domain, and the cross-correlation function data segments between each station pair were calculated;

[0006] The root mean square ratio screening and superposition method is used to screen the cross-correlation function data segments;

[0007] The selected multichannel cross-correlation function data segments are transformed to the frequency domain to construct the cross-spectrum matrix at each frequency point;

[0008] Based on the cross-spectral matrix, design an adaptive filter;

[0009] The cross-correlation function data segment is filtered using the adaptive filter to obtain the filtered signal;

[0010] The filtered signal is decomposed into time and frequency components, and the instantaneous phase is extracted.

[0011] A time-frequency domain coherent weight is constructed based on the instantaneous phase information;

[0012] The filtered signal is weighted and superimposed using the time-frequency domain coherence weights to obtain the superimposed time-frequency spectrum.

[0013] The superimposed time-frequency spectrum is subjected to inverse time-frequency transformation to obtain the final time-domain superposition result.

[0014] Furthermore, the cross-correlation function data segments between each station pair are calculated in the frequency domain, including:

[0015] Acquire the background noise time-domain signal synchronously recorded by two stations, one of which is mathematically equivalent to a virtual source and the other station is a receiver;

[0016] The time-domain signals from the two stations are transformed to the frequency domain using Fourier transform to obtain the corresponding spectra.

[0017] In the frequency domain, the spectrum of one station is taken as its complex conjugate and multiplied with the spectrum of the other station to obtain the frequency domain cross-spectrum result.

[0018] Perform an inverse Fourier transform on the frequency domain cross-spectral results and transform them back to the time domain to obtain the cross-correlation function data segment between the two stations.

[0019] Furthermore, the root mean square ratio screening and superposition method is used to screen the cross-correlation function data segments, including:

[0020] Calculate the linear superposition result of all cross-correlation function data segments, and use its signal-to-noise ratio as the benchmark reference value;

[0021] For the Given a cross-correlation function data segment, construct a temporary linear superposition result of the remaining cross-correlation function data segments excluding the given cross-correlation function data segment;

[0022] Calculate the root mean square (RMS) values ​​of the temporary linear superposition result in the positive delay signal window, the negative delay signal window, and the noise window with positive delay and negative delay, respectively.

[0023] Based on the root mean square value, calculate the value excluding the first... The second signal-to-noise ratio of the remaining cross-correlation function data segments outside of the first cross-correlation function data segments;

[0024] The second signal-to-noise ratio and the benchmark reference value are used as quality evaluation factors.

[0025] Based on the comparison results between the quality assessment factor and the dynamic threshold, the required set of cross-correlation function data segments is selected.

[0026] Furthermore, based on the comparison results between the quality assessment factor and the dynamic threshold, the required set of cross-correlation function data segments is selected, including:

[0027] Obtain the total number of cross-correlation function data segments ;

[0028] Based on the total number Determine the dynamic threshold The dynamic threshold With the total quantity Positive correlation, such that when the total number When it increases, the dynamic threshold The dynamic threshold T approaches 1 to implement stricter screening, and decreases accordingly to implement more lenient screening as the total number N decreases; and the dynamic threshold T is used to screen multiple cross-correlation function data segments.

[0029] Furthermore, the root mean square ratio screening and superposition method is used to screen the cross-correlation function data segments, which also includes:

[0030] Based on the group velocity model, the positive and negative time delay intervals corresponding to the effective seismic wave signal are determined. The positive time delay interval corresponds to the time delay range of the seismic wave propagating from the virtual source to the receiver in the forward direction, and the negative time delay interval corresponds to the time delay range of the seismic wave propagating from the receiver to the virtual source in the reverse direction.

[0031] Within the positive time delay interval, a continuous time interval covering the time domain range of the effective seismic wave signal is selected as the positive time delay signal window; within the negative time delay interval, a continuous time interval covering the time domain range of the effective seismic wave signal is selected as the negative time delay signal window.

[0032] Within the positive delay interval, two consecutive time intervals without valid signals are selected on both sides of the positive delay signal window, respectively, as the first noise window and the second noise window corresponding to the positive delay; within the negative delay interval, two consecutive time intervals without valid signals are selected on both sides of the negative delay signal window, respectively, as the first noise window and the second noise window corresponding to the negative delay.

[0033] Furthermore, the cross-spectral matrix refers to a two-dimensional matrix formed by arranging the cross-spectral results of all cross-correlation function pairs at frequency points according to the rule of pairing corresponding signals in rows and columns, using the filtered cross-correlation function data segments as indices.

[0034] Furthermore, the filter coefficients of the adaptive filter Represented as:

[0035] ,

[0036] In the formula, This indicates the number of cross-correlation function data segments retained after filtering. This indicates that all signal channels are paired at the frequency point. The corresponding cross spectrum matrix element Perform full element summation.

[0037] This indicates that all diagonal elements are at frequency points. The values ​​at each location are accumulated. Indicates the first One signal channel, Indicates the first One signal channel.

[0038] Furthermore, the filtered signal The calculation formula is:

[0039] ,

[0040] In the formula, The parameter used to adjust the roughness of the filter has a value greater than 0. This is a data segment representing the cross-correlation function in the frequency domain. To delay the transmission time, , For the index of stations, For time.

[0041] A seismic interferometric signal enhancement system based on an adaptive superposition strategy is provided in the second aspect of this application, comprising:

[0042] The cross-correlation calculation module is used to convert the seismic interferometric time-domain signal to the frequency domain and calculate the cross-correlation function data segment;

[0043] The quality screening module is used to screen cross-correlation function data segments using the root mean square ratio screening superposition method;

[0044] The cross-spectral analysis module is used to construct the cross-spectral matrix of the filtered data and design adaptive filters;

[0045] The filtering module is used to filter the cross-correlation function data segment using an adaptive filter to obtain the filtered signal.

[0046] The time-frequency analysis module is used to perform time-frequency decomposition on the filtered signal and extract instantaneous phase information in the frequency domain.

[0047] The weight construction module is used to construct time-frequency domain coherent weights based on instantaneous phase information.

[0048] The weighted superposition module is used to weight and superimpose the filtered signal using time-frequency domain coherent weights to obtain the superimposed time-frequency spectrum.

[0049] The inverse transform module is used to inversely transform the superimposed time spectrum back to the time domain to obtain the final time-domain superposition result.

[0050] The embodiments of this application have at least the following beneficial effects:

[0051] This application can effectively reduce noise interference in signals, making the signal waveform more stable and regular. Compared with existing single / basic signal processing methods, it can achieve noise reduction and enhancement of target signals more efficiently, significantly improving signal quality. Attached Figure Description

[0052] Figure 1 This is a flowchart of the seismic interferometric signal enhancement method based on an adaptive stacking strategy provided in the embodiments of this application;

[0053] Figure 2 This is the spatial distribution of stations and noise sources provided in the embodiments of this application;

[0054] Figure 3 The results of superimposing cross-correlation functions are obtained by different superposition methods provided in the embodiments of this application;

[0055] Figure 4 This is a spatial distribution map of the sources retained and eliminated by the method provided in the embodiments of this application. (a) is the spatial distribution of the retained sources, and (b) is the spatial distribution of the eliminated sources.

[0056] Figure 5 This is the original seismic waveform record diagram provided in the embodiments of this application;

[0057] Figure 6 These are the dispersion energy spectra of the station pair cross-correlation function provided in the embodiments of this application; wherein, (a) is the group velocity dispersion spectrum of the linear superposition method, (b) is the group velocity dispersion spectrum of the autocorrelation function + phase-weighted superposition method, (c) is the group velocity dispersion spectrum of the autocorrelation function + time-frequency domain phase-weighted superposition method, (d) is the phase velocity dispersion spectrum of the linear superposition method, (e) is the phase velocity dispersion spectrum of the autocorrelation function + phase-weighted superposition method, and (f) is the phase velocity dispersion spectrum of the autocorrelation function + time-frequency domain phase-weighted superposition method. (g) is the group velocity dispersion diagram of the phase velocity dispersion diagram of the phase velocity dispersion diagram of the root mean square ratio screening method, (h) is the velocity dispersion diagram of the autocorrelation function based on the root mean square ratio + phase velocity dispersion diagram of the phase velocity dispersion diagram of the phase velocity dispersion diagram of the method of this application, (j) is the phase velocity dispersion diagram of the root mean square ratio screening method, (k) is the phase velocity dispersion diagram of the autocorrelation function based on the root mean square ratio + phase velocity dispersion diagram of the phase velocity dispersion diagram of the method of this application, and (l) is the phase velocity dispersion diagram of the method of this application.

[0058] Figure 7 The actual measured average dispersion data provided in the embodiments of this application are compared with the theoretical dispersion curve obtained by forward modeling using an initial one-dimensional velocity model;

[0059] Figure 8 This is a horizontal cross-sectional view of the S-wave velocity structure at a depth of 0.6 km provided in an embodiment of this application;

[0060] Figure 9 This is a seismic reflection profile provided in an embodiment of this application;

[0061] Figure 10 This is a horizontal cross-sectional view of the S-wave velocity structure with the earthquake reflection profile superimposed, provided in an embodiment of this application. Detailed Implementation

[0062] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0063] The seismic interferometric signal enhancement method and system based on the adaptive stacking strategy provided in this application firstly identify cross-correlation functions (CCFs) that make constructive contributions to the reconstruction of the empirical Green's function (EGF) based on the screening criterion of the root mean square ratio screening stacking method, thus solving the problem of coherent signal identification; secondly, by using adaptive covariance filtering and time-frequency domain phase weighted stacking techniques, the instantaneous phase coherence is improved in multiple frequency bandwidths while suppressing random phase noise components.

[0064] See Figure 1 As shown, a seismic interferometric signal enhancement method based on an adaptive stacking strategy includes:

[0065] S101, convert the seismic interferometric time-domain signal to the frequency domain and calculate the cross-correlation function data segment between each station pair;

[0066] S102, the root mean square ratio screening and superposition method is used to screen the cross-correlation function data segments;

[0067] S103, transform the selected multi-channel cross-correlation function data segments to the frequency domain and construct the cross-spectrum matrix at each frequency point;

[0068] S104, Design an adaptive filter based on the cross-spectral matrix;

[0069] S105, The cross-correlation function data segment is filtered using the adaptive filter to obtain the filtered signal;

[0070] S106, Perform time-frequency decomposition on the filtered signal and extract instantaneous phase information;

[0071] S107, Construct time-frequency domain coherent weights based on the instantaneous phase information;

[0072] S108, The filtered signal is weighted and superimposed using the time-frequency domain coherent weights to obtain the superimposed time-frequency spectrum;

[0073] S109, Perform inverse time-frequency transformation on the superimposed time spectrum to obtain the final time-domain superposition result.

[0074] In step S101, the seismic interferometric time-domain signal is converted to the frequency domain using Fourier transform. The cross-correlation function data segments between each station pair are then calculated. Frequency domain cross-correlation is the core computational method for obtaining the propagation delay between the virtual source and receiver. The peak position of the cross-correlation function data segment corresponds to the propagation delay. In the frequency domain cross-correlation operation, the multiplication process of the virtual source signal's frequency domain complex conjugate and the receiver signal's frequency domain automatically extracts the phase difference between the two. , , These are the frequency domain phases of the virtual source signal and the receiver signal, respectively; this phase difference is related to the propagation delay from the virtual source to the receiver. satisfy ( (where angular frequency is the metric), and the peak position of the cross-correlation function data segment of the frequency domain cross-correlation output corresponds to this propagation delay. The data segment for calculating the cross-correlation function between each station pair in the frequency domain includes:

[0075] Acquire the background noise time-domain signal synchronously recorded by two stations, one of which is mathematically equivalent to a virtual source and the other station is a receiver;

[0076] The time-domain signals from the two stations are transformed to the frequency domain using Fourier transform to obtain the corresponding spectra.

[0077] In the frequency domain, the spectrum of one station is taken as its complex conjugate and multiplied with the spectrum of the other station to obtain the frequency domain cross-spectrum result.

[0078] Perform an inverse Fourier transform on the frequency domain cross-spectral results and transform them back to the time domain to obtain the cross-correlation function data segment between the two stations.

[0079] The cross-correlation function data segment of the frequency is calculated using the following formula:

[0080] ,

[0081] In the formula, This represents the cross-correlation function data segment between the virtual source and receiver, where... and These are the background noise time-domain signals recorded by two stations. Although physically it is a receiving point, its background noise time-domain signal is mathematically equivalent to a virtual seismic source, while the station... It still acts as a receiver. This represents the Fourier transform operator. This represents the complex conjugate of the Fourier transform. This represents the inverse Fourier transform operator.

[0082] In step S102, when the coherent signal is overwhelmed by strong random noise, superimposing only high-quality cross-correlation function data segments and removing low-quality cross-correlation function data segments helps to recover the Green's function. Based on the group velocity model, the group velocity to time interval of the path between receiver pairs in the period of interest is calculated. A signal window with positive time delay is defined, and two noise windows are defined respectively. Similarly, a signal window and two noise windows are symmetrically set for the negative time delay of the data.

[0083] Specifically, based on the swarm velocity model, the positive and negative time delay intervals corresponding to the effective seismic wave signal are determined. The positive time delay interval corresponds to the time delay range of the seismic wave propagating from the virtual source to the receiver in the forward direction, while the negative time delay interval corresponds to the time delay range of the seismic wave propagating from the receiver to the virtual source in the reverse direction. The swarm velocity model is prior knowledge, such as previous exploration, sonic logging, or velocity models from previous researchers in the same study area. It is a mathematical model describing the variation of the seismic wave swarm velocity (i.e., the speed of seismic wave energy propagation) with parameters such as propagation medium properties, propagation distance, and seismic wave frequency. This swarm velocity model is based on the propagation mechanism of seismic waves in the underground medium and is obtained by fitting geological exploration experience data with measured seismic wave data. Its core function is to determine the propagation time interval of the seismic wave signal in the time domain.

[0084] Within the positive time delay interval, a continuous time interval covering the effective time domain range of the seismic wave signal is selected as the positive time delay signal window; similarly, within the negative time delay interval, a continuous time interval covering the effective time domain range of the seismic wave signal is selected as the negative time delay signal window. The boundary determination of the signal window must comprehensively consider the dominant frequency of the seismic wave, the medium quality factor, and the geometric diffusion effect of the expected propagation path to ensure that most of the effective signal energy is completely captured. For the positive time delay signal window, its starting time is typically set to the theoretical arrival time of the first seismic wave after virtual source excitation minus a small tolerance to compensate for the uncertainty of the group velocity model; the ending time is determined based on the expected arrival time of the farthest effective reflected or scattered signal, while avoiding the inclusion of subsequent strongly coherent noise into the window. The setting of the negative time delay signal window follows the principle of complete symmetry, that is, with zero time delay as the center, the positive time delay signal window is mirrored and flipped along the time axis to the negative time delay interval, thereby ensuring the statistical consistency of the positive and negative time delay analysis framework.

[0085] Within the positive delay interval, two consecutive time intervals without valid signals are selected on both sides of the positive delay signal window, serving as the first noise window and the second noise window corresponding to the positive delay; within the negative delay interval, two consecutive time intervals without valid signals are selected on both sides of the negative delay signal window, serving as the first noise window and the second noise window corresponding to the negative delay.

[0086] In the regions immediately adjacent to the positive and negative time-delay signal windows, a first noise window and a second noise window are respectively set. The first noise window is located before the positive time-delay signal window and is used to characterize the preceding background noise and potential incoherent interference; the second noise window is located after the positive time-delay signal window and is used to capture trailing noise and possible residual energy of multiple waves. The widths of the two noise windows are usually comparable to or slightly narrower than the corresponding signal windows, and an appropriate time gap is maintained between them to avoid signal leakage contaminating the noise statistical characteristics. Specifically, the first noise window of the positive time-delay signal window within the positive time-delay interval starts at the starting boundary of the positive time-delay interval and ends at the starting time of the positive time-delay signal window minus the guard gap; the second noise window starts at the ending time of the positive time-delay signal window plus the guard gap and ends at the ending boundary of the positive time-delay interval. The noise windows on both sides of the negative time-delay signal window within the negative time-delay interval are set strictly symmetrically with those of the positive time-delay interval. The first noise window corresponding to the negative time-delay corresponds to the preceding noise propagating in the backward direction, and the second noise window corresponds to the trailing noise propagating in the backward direction.

[0087] The core purpose of this windowing strategy is to establish a quantitative separation framework for signal and noise, providing a reliable statistical basis for subsequent adaptive superposition weight calculation. By comparing the energy distribution characteristics within the signal window and the noise window, it is possible to effectively identify and suppress recording channels with low signal-to-noise ratio and poor coherence, while enhancing the contribution weight of high signal-to-noise ratio channels.

[0088] The purpose of the above is to filter cross-correlation function data segments using the root mean square ratio screening and superposition method, specifically including:

[0089] Calculate the linear superposition result of all cross-correlation function data segments, and use its signal-to-noise ratio as the benchmark reference value;

[0090] For the Given a cross-correlation function data segment, construct a temporary linear superposition result of the remaining cross-correlation function data segments excluding the given cross-correlation function data segment;

[0091] Calculate the root mean square (RMS) values ​​of the temporary linear superposition result in the positive delay signal window, the negative delay signal window, and the noise window with positive delay and negative delay, respectively.

[0092] Based on the root mean square value, calculate the value excluding the first... The second signal-to-noise ratio of the remaining cross-correlation function data segments outside of the first cross-correlation function data segments;

[0093] The second signal-to-noise ratio and the benchmark reference value are used as quality evaluation factors.

[0094] Based on the comparison results between the quality assessment factor and the dynamic threshold, the required set of cross-correlation function data segments is selected.

[0095] By screening, cross-correlation function data segments that contribute positively to the recovery of the Green's function are identified.

[0096] Calculate the root mean square (RMS) value of the window. Calculate the RMS value of each cross-correlation function data segment within the positive time-delay signal window, the negative time-delay signal window, and each noise window. The formula is as follows:

[0097] ,

[0098] ,

[0099] In the formula, and These represent the positive delay signal window and the negative delay signal window, respectively. and These represent two noise windows with positive time delay and two noise windows with negative time delay, respectively. Represents the root mean square function. Represents the root mean square value of the signal. This represents the root mean square value of the noise.

[0100] The signal-to-noise ratio of the linear superposition of all cross-correlation function data segments is used as a benchmark:

[0101] ,

[0102] in, This is represented as a baseline reference value.

[0103] For each cross-correlation function data segment, a temporary linear superposition result is calculated. The temporary linear superposition result is excluding the first segment. The superposition of the remaining cross-correlation function data segments after the first cross-correlation function data segment, where superposition refers to the superposition of the remaining cross-correlation function data segments after the first cross-correlation function data segment. Except for the first cross-correlation function data segment, all remaining cross-correlation function data segments, after pre-signal processing, noise window division, and preliminary quality screening, are subjected to linear superposition in the time domain. The superposition process preserves the time-domain characteristics (such as phase and amplitude trends) of each data segment; similarly, calculations are performed except for the first cross-correlation function data segment. Signal-to-noise ratio after each cross-correlation function data segment Here, the second signal-to-noise ratio is used, and the calculation method for the second signal-to-noise ratio is: based on the first... In addition to the cross-correlation function data segments, the root mean square values ​​(denoted as ) of the positive and negative time delay signal windows corresponding to all remaining cross-correlation function data segments are also included. ) and the root mean square values ​​of the noise windows corresponding to positive and negative time delays (denoted as ) According to the formula Calculated; where, This is the arithmetic mean of the root mean square values ​​of the remaining segment within the positive and negative time delay signal windows. The arithmetic mean of the root mean square values ​​of the cross-correlation function data segment within the corresponding positive and negative time delay noise windows; the second signal-to-noise ratio and the benchmark reference value are defined as quality assessment factors:

[0104] ,

[0105] Therefore, considering both positive and negative time delays, the decision logic is designed as follows:

[0106] ,

[0107] ,

[0108] In the formula, Represents a dynamic threshold. Representing the number of cross-correlation function segments, based on the comparison results of the quality assessment factor and the dynamic threshold, the required set of cross-correlation function data segments is selected, including:

[0109] Obtain the total number N of cross-correlation function data segments;

[0110] Based on the total number Determine the dynamic threshold The dynamic threshold With the total quantity Positive correlation, such that when the total number When it increases, the dynamic threshold Approaching 1 to implement more stringent screening, when the total number When it decreases, the dynamic threshold Accordingly, the threshold is lowered to implement a more lenient screening; and the dynamic threshold is utilized. We can filter multiple cross-correlation function data segments to obtain a high-quality set of cross-correlation functions.

[0111] Steps 103 to 105 involve designing an adaptive filter for the filtered multi-channel cross-correlation function data segment and filtering it. Specifically, this includes:

[0112] After obtaining a high-quality set of cross-correlation function data segments, these filtered cross-correlation function data segments are processed. Noise is further suppressed by enhancing the coherent signals on all filtered cross-correlation function data segments. Specifically, for each pair of filtered cross-correlation function data segments, they are converted to a frequency domain representation;

[0113] The cross-spectral matrix is ​​calculated as follows:

[0114] ,

[0115] in, For the frequency domain The cross-correlation function data segment of each signal channel, For the frequency domain The cross-correlation function data segment of each signal channel. for One signal channel and The cross-spectral matrix of the cross-correlation function data segment of each signal channel is equivalent to the time covariance matrix.

[0116] The cross-spectral matrix mentioned here refers to a two-dimensional matrix formed by arranging the cross-spectral results of all cross-correlation function data segments at a given frequency point according to the rule of "row and column corresponding signal pairing" using the pre-selected cross-correlation function data segments as indices. For example, if there are A pre-selected cross-correlation function data segments, the cross-spectral matrix is ​​A×A dimensional. The element in the a-th row and b-th column of the cross-spectral matrix corresponds to the cross-spectral results of the a-th and b-th cross-correlation function data segments at the current frequency point, reflecting the amplitude and phase correlation relationship between the corresponding two cross-correlation function data segments at that frequency point. This cross-spectral matrix reflects the mutual relationship of each cross-correlation function data segment at that frequency point. Based on the constructed cross-spectral matrix, a filter coefficient is designed for each pair of cross-correlation function data segments. The design principle of the filter coefficient is: when the two cross-correlation function data segments are highly coherent, the filter coefficient approaches 1, allowing the signal to pass through with almost no loss; when the two cross-correlation function data segments are uncorrelated, the coefficient approaches 0, thereby effectively suppressing noise components. Therefore, the filter can adaptively adjust between 0 and 1 according to the degree of coherence between the cross-correlation function data segments.

[0117] The filter coefficients of the adaptive filter at each frequency Represented as:

[0118] ,

[0119] In the formula, This indicates the number of cross-correlation function data segments retained after filtering. This indicates that all signal channels are paired at the frequency point. The corresponding cross spectrum matrix element Perform full element summation.

[0120] This indicates that all diagonal elements are at frequency points. The values ​​at each location are accumulated. Indicates the first One signal channel, Indicates the first One signal channel.

[0121] Filtered signal The calculation formula is:

[0122] ,

[0123] In the formula, The parameter used to adjust the roughness of the filter has a value greater than 0. This is a data segment representing the cross-correlation function in the frequency domain. To delay the transmission time, , For the index of stations, For time. By introducing an adjustable power parameter. This controls the "roughness" or sharpness of the filter, thereby controlling the noise reduction intensity.

[0124] Steps 106 to 109 involve performing time-frequency decomposition on the filtered signal and extracting the instantaneous phase; constructing time-frequency domain coherent weights based on the instantaneous phase information; using the time-frequency domain coherent weights to perform weighted superposition on the filtered signal to obtain the superimposed time spectrum; and performing inverse time-frequency transformation on the superimposed time spectrum to obtain the final time-domain superposition result.

[0125] The time-frequency decomposition of each filtered signal is obtained based on the S-transform:

[0126] ,

[0127] in, Therefore A Gaussian window function centered on the center, with width and Proportional For frequency, The propagation delay is the position corresponding to the peak position of the cross-correlation function data segment. This is the filtered signal. The time-frequency decomposition signal is the time-frequency decomposition signal of each filtered signal, which can be represented by instantaneous amplitude and instantaneous phase.

[0128] The time-frequency domain coherence weights are obtained by defining the instantaneous phase sum:

[0129] ,

[0130] in, For time-frequency domain coherence weights, the exponent Control the suppression intensity of incoherent components.

[0131] The filtered signals are weighted and superimposed to obtain the superimposed time spectrum, which is expressed as:

[0132] ;

[0133] Performing an inverse S-transform on the above formula yields the final time-domain superposition result:

[0134] .

[0135] On the other hand, embodiments of this application provide a seismic interferometric signal enhancement system based on an adaptive stacking strategy, comprising:

[0136] The cross-correlation calculation module is used to convert the seismic interferometric time-domain signal to the frequency domain and calculate the cross-correlation function data segment;

[0137] The quality screening module is used to screen cross-correlation function data segments using the root mean square ratio screening superposition method;

[0138] The cross-spectral analysis module is used to construct the cross-spectral matrix of the filtered data and design adaptive filters;

[0139] The filtering module is used to filter the cross-correlation function data segment using an adaptive filter to obtain the filtered signal.

[0140] The time-frequency analysis module is used to perform time-frequency decomposition on the filtered signal and extract instantaneous phase information in the frequency domain.

[0141] The weight construction module is used to construct time-frequency domain coherent weights based on instantaneous phase information.

[0142] The weighted superposition module is used to weight and superimpose the filtered signal using time-frequency domain coherent weights to obtain the superimposed time-frequency spectrum.

[0143] The inverse transform module is used to inversely transform the superimposed time spectrum back to the time domain to obtain the final time-domain superposition result.

[0144] To quantitatively evaluate the effectiveness of the proposed method in suppressing high-density noise sources in the unstable phase region, the following numerical experiments were designed:

[0145] In the experimental setup, two centrally symmetrical stations were deployed, with a spacing of 150 km between them. Within a ring-shaped region 200-600 km from the center of the array, approximately 18,000 noise sources were randomly distributed within a vertical depth of 0-30 km to simulate the regional background noise field. To simulate the influence of highly directional noise sources, concentrated noise sources with three times the density were arranged within a fan-shaped region at an azimuth angle of 45°-60° and a radial distance of 500-600 km, with their spatial distribution as shown below. Figure 2 As shown in (a), Figure 2 (b) in the figure is a waveform recording diagram of the signals from station A and station B. The vertical axis represents the amplitude, which is used to show the dynamic changes of the signals collected by different stations over time.

[0146] All noise sources employed Ricker wavelets with a dominant frequency range of 0.05-0.2 Hz and a uniform amplitude distribution between 0.5 and 1 Hz, superimposed with Gaussian white noise of mean 0 and standard deviation 0.05. A dynamic noise source distribution model was used: starting from a 30 azimuth angle, the source rotated counterclockwise by 5° every 10,000 s, with a complete rotation period of 720,000 s, divided into 72 time intervals. The precise coordinates and excitation time of each noise source were recorded in detail during the experiment to facilitate subsequent quantitative analysis of the source location's impact on the Green's function reconstruction quality.

[0147] The subsurface medium model adopts a three-layer horizontal stratified model based on a simplified three-dimensional S-wave velocity structure. According to empirical formulas, the corresponding P-wave velocity and density parameters are derived from the S-wave velocity, and the specific parameters are shown in Table 1.

[0148] Table 1. Parameters of the three-layer velocity model:

[0149] Floor number Thickness (km) Shear wave velocity (km / s) Longitudinal wave velocity (km / s) Density (g / m^3) 1 10 2.89 4.88 2.52 2 20 3.47 5.91 2.70 3 Infinity 3.74 6.43 2.83

[0150] First, the original seismic waveform data is segmented into 1000-second segments. After data processing, frequency domain cross-correlation is calculated on the segmented waveforms. To evaluate the effectiveness of different stacking methods, the proposed method is compared with other stacking methods. The cross-correlation function after stacking is shown below. Figure 3 As shown. Figure 3 The figure shows a comparison of the signal time-domain waveforms and signal-to-noise ratio (SNR) performance of different signal processing methods. This figure can intuitively demonstrate the signal processing effect of the technical solution of this application. Specifically, the vertical axis corresponds to six signal processing methods, including the linear processing method, the RMSR method (root mean square ratio), the ACF+PWS combination method (autocorrelation function + phase-weighted superposition combination method), the ACF+tf-PWS combination method (autocorrelation function + time-frequency domain phase-weighted superposition combination method), the RMSR+ACF+PWS combination method (root mean square ratio + autocorrelation function + time-frequency domain phase-weighted superposition combination method, the method of this application); the horizontal axis represents time (unit: s). Figure 3 The curves in the middle represent the time-domain waveforms of the signals corresponding to each processing method; the "SNR-neg" (negative signal-to-noise ratio) and "SNR-pos" (positive signal-to-noise ratio) labels next to each processing method are quantitative indicators of signal quality. Higher values ​​of these indicators indicate better noise reduction and enhancement effects. Figure 3Quantitative data and waveform characteristics show that, compared to a single linear processing method (with an SNR-neg of only 9.6 and an SNR-pos of 8.5), the combined processing method proposed in this application significantly improves the signal-to-noise ratio, achieving an SNR-neg of 28.1 and an SNR-pos of 27.5. Furthermore, considering the fluctuation regularity of the time-domain waveform, it is evident that the combined processing method of this application effectively reduces noise interference in the signal, making the signal waveform more stable and regular. Therefore, it is clear that this application can more efficiently achieve noise reduction and enhancement of the target signal, significantly improving signal quality and thus ensuring the accuracy of subsequent applications based on this signal.

[0151] The signal enhancement effect of each method is quantitatively evaluated by calculating the signal-to-noise ratio (SNR). The calculation formula is as follows:

[0152] ,

[0153] In the formula, It is the root mean square operator; and The amplitudes of the waveforms within the signal window and trailing noise window are respectively. The signal-to-noise ratio (SNR) is calculated and analyzed for the positive and negative delay branches of the cross-correlation function data segment. Experimental results show that different superposition strategies have significant differences in improving the quality of the empirical Green's function. The traditional linear superposition method only achieves an average SNR of 8.95 dB, while the method proposed in this application achieves an optimal SNR of 26.9 dB.

[0154] Figure 4 (a) and Figure 4 Figure (b) shows the spatial distribution of the retained and removed sources, respectively. The results indicate that the method described in this application can effectively distinguish between sources that contribute positively to the recovery of the empirical Green's function and sources that introduce interference. Figure 4 In diagram (b), the red hyperbolic dashed line represents the boundary of the stable phase region. The formula for calculating the angle of the stable phase region is:

[0155] ,

[0156] in It's the wavelength. It is the length of the station's continuous line. For the stable phase region angle. Figure 4 In (a) of the model, the stable-phase region along the receiver path exhibits a significantly higher source density, and most sources located within this region are preserved during superposition. These stable-phase regions are crucial for the coherent recovery of the empirical Green's function. Conversely, Figure 4(b) shows that most of the rejected sources clustered in the unstable phase region, which introduces spurious arrivals and degrades signal quality in the cross-correlation function data segment. Notably, most sources outside the stable phase region were rejected, while those within the stable phase region were retained. Even in regions with increased source density, the method of this application maintains its discriminative capability. Here, SPZ refers to potential source zones, which are typically the focus of medium- to long-term seismic hazard assessments.

[0157] To verify the effectiveness and accuracy of the method proposed in this application, a field observation experiment was conducted in an area with complete geological background data. The experimental area was located in a certain region, and the observation period was from October 30 to November 20, 2024. A dense array of 34 three-component seismographs was deployed. The instruments were broadband seismographs, and the array coverage area was approximately 6 km × 7 km. The data acquisition sampling rate was 1000 Hz.

[0158] Figure 5 The original seismic waveform record is displayed. Figure 5 The stations in the selected area received seismic waves during the corresponding time period, which are the valid signal segments corresponding to the seismic events captured in this observation.

[0159] In the data preprocessing stage, the raw data is first downsampled to 50Hz. Cross-correlation in the frequency domain is used for calculation, and the waveforms of cross-correlation functions obtained by different superposition methods are compared. Compared with the traditional linear superposition method, this application shows a significant advantage in suppressing incoherent noise.

[0160] To further verify the effectiveness of the method, stations 786-809 were selected for analysis. Figure 6 The dispersion energy spectrum of the cross-correlation function of the stations is given, where, Figure 6 (a), (b), and (c) in the text. Figure 6 (g), (h), and (i) in the diagram represent the group velocity dispersion plot; Figure 6 (d), (e), and (f) in the text. Figure 6 In this context, (j), (k), and (l) represent phase velocity dispersion diagrams. Different superposition methods include: Figure 6 In the diagram, (a) and (d) are examples of linear superposition. Figure 6 In the middle (b) and (e), the autocorrelation function is used in the phase-weighted superposition method; Figure 6 The autocorrelation functions (c) and (f) in the middle are combined with the time-frequency domain phase-weighted superposition method; Figure 6 (g) and (j) in the figure represent the root mean square ratio screening method; Figure 6 In this context, (h) and (k) are autocorrelation functions based on the root mean square ratio plus phase-weighted superposition method; Figure 6In the diagram, (i) and (l) are based on the method of this application, and the black and white dotted lines represent the extracted group velocity and phase velocity, respectively. The comparative results show that the incoherent noise component in the cross-correlation function data after processing by the method of this application is significantly reduced, and the dispersion energy spectrum exhibits higher frequency resolution, providing a high-quality data foundation for subsequent surface wave dispersion measurements.

[0161] In the dispersion curve extraction stage, a convolutional neural network method is used to automatically extract phase velocity information from the processed spectral data. The semi-analytical spectral element method is then used to perform forward modeling of the theoretical dispersion curve of the initial one-dimensional velocity model. Figure 7 The results show that the theoretical dispersion curve calculated by the theoretical model has good fitting accuracy with the observed average dispersion data. Finally, a three-dimensional velocity structure model of the S-wave in the study area was constructed using a direct inversion method. Figure 8 A horizontal profile of the S-wave velocity structure at a depth of 0.6 km is presented. To verify the accuracy and reliability of the inverted S-wave velocity structure, it is compared and analyzed with seismic reflection profile data of the study area. Figure 9 This is a seismic reflection profile, with solid green lines indicating the inferred fault locations. Figure 10 The S-wave velocity structure superimposed on the seismic reflection profile demonstrates the superposition results of the vertical profile of the S-wave velocity model and the seismic reflection record. Comparative analysis shows that the spatial distribution of the inverted S-wave velocity structure is in good agreement with the subsurface structural features shown by the seismic reflection profile, confirming the effectiveness of the method proposed in this application and thus ensuring the geological rationality and reliability of the tomographic imaging results.

[0162] Based on the above comparisons and field verifications, it is evident that the method presented in this application effectively eliminates unstable phase interference, significantly improves the signal-to-noise ratio (SNR), and practical applications demonstrate that this method requires only about ten days of noise data to obtain cross-correlation function data segments with high SNR, resulting in a more concentrated and continuous dispersion energy distribution and greatly improving the extraction quality and computational efficiency of surface wave dispersion curves. In summary, the method presented in this application demonstrates excellent performance in improving the SNR of cross-correlation functions and the quality of dispersion analysis, and has broad application prospects in seismic exploration, subsurface structure imaging, and other fields, particularly suitable for high-precision seismic analysis under short-term observation data conditions.

[0163] The above description is merely a preferred embodiment of this application and is not intended to limit this application. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this application should be included within the protection scope of this application.

Claims

1. A seismic interferometric signal enhancement method based on an adaptive stacking strategy, characterized in that, The method includes: The seismic interferometric time-domain signal was converted to the frequency domain, and the cross-correlation function data segments between each station pair were calculated; The root mean square ratio screening and superposition method is used to screen the cross-correlation function data segments; The selected multichannel cross-correlation function data segments are transformed to the frequency domain to construct the cross-spectrum matrix at each frequency point; Based on the cross-spectral matrix, design an adaptive filter; the filter coefficients of the adaptive filter... Represented as: ， In the formula, This indicates the number of cross-correlation function data segments retained after filtering. This indicates that all signal channels are paired at the frequency point. The corresponding cross spectrum matrix element Perform full element-wise summation. This indicates that all diagonal elements are at frequency points. The values ​​at each location are accumulated. Indicates the first One signal channel, Indicates the first One signal channel; The cross-correlation function data segment is filtered using the adaptive filter to obtain the filtered signal; The filtered signal is decomposed into time and frequency components, and the instantaneous phase is extracted. Construct time-frequency domain coherent weights based on instantaneous phase information; The filtered signal is weighted and superimposed using the time-frequency domain coherence weights to obtain the superimposed time-frequency spectrum. The superimposed time-frequency spectrum is subjected to inverse time-frequency transformation to obtain the final time-domain superposition result.

2. The seismic interferometric signal enhancement method based on an adaptive superposition strategy according to claim 1, characterized in that, The data segment used to calculate the cross-correlation function between each station pair in the frequency domain includes: Acquire the background noise time-domain signal synchronously recorded by two stations, one of which is mathematically equivalent to a virtual source and the other station is a receiver; The time-domain signals from the two stations are transformed to the frequency domain using Fourier transform to obtain the corresponding spectra. In the frequency domain, the spectrum of one station is taken as its complex conjugate and multiplied with the spectrum of the other station to obtain the frequency domain cross-spectrum result. Perform an inverse Fourier transform on the frequency domain cross-spectral results and transform them back to the time domain to obtain the cross-correlation function data segment between the two stations.

3. The seismic interferometric signal enhancement method based on an adaptive stacking strategy according to claim 2, characterized in that, The root mean square ratio screening and superposition method is used to screen cross-correlation function data segments, including: Calculate the linear superposition result of all cross-correlation function data segments, and use its signal-to-noise ratio as the benchmark reference value; For the Given a cross-correlation function data segment, construct a temporary linear superposition result of the remaining cross-correlation function data segments excluding the given cross-correlation function data segment; Calculate the root mean square (RMS) values ​​of the temporary linear superposition result in the positive delay signal window, the negative delay signal window, and the noise window with positive delay and negative delay, respectively. Based on the root mean square value, calculate the value excluding the first... The second signal-to-noise ratio of the remaining cross-correlation function data segments outside of the first cross-correlation function data segments; The second signal-to-noise ratio and the benchmark reference value are used as quality evaluation factors. Based on the comparison results between the quality assessment factor and the dynamic threshold, the required set of cross-correlation function data segments is selected.

4. The seismic interferometric signal enhancement method based on an adaptive superposition strategy according to claim 3, characterized in that, Based on the comparison results between the quality assessment factor and the dynamic threshold, the required set of cross-correlation function data segments is selected, including: Obtain the total number of cross-correlation function data segments ; Based on the total number Determine the dynamic threshold The dynamic threshold With the total quantity Positive correlation, such that when the total number When it increases, the dynamic threshold The dynamic threshold T approaches 1 to implement stricter screening, and decreases accordingly to implement more lenient screening as the total number N decreases; and the dynamic threshold T is used to screen multiple cross-correlation function data segments.

5. The seismic interferometric signal enhancement method based on an adaptive superposition strategy according to claim 3, characterized in that, The root mean square ratio screening and superposition method is used to screen cross-correlation function data segments, and also includes: Based on the group velocity model, the positive and negative time delay intervals corresponding to the effective seismic wave signal are determined. The positive time delay interval corresponds to the time delay range of the seismic wave propagating from the virtual source to the receiver in the forward direction, and the negative time delay interval corresponds to the time delay range of the seismic wave propagating from the receiver to the virtual source in the reverse direction. Within the positive time delay interval, a continuous time interval covering the time domain range of the effective seismic wave signal is selected as the positive time delay signal window; within the negative time delay interval, a continuous time interval covering the time domain range of the effective seismic wave signal is selected as the negative time delay signal window. Within the positive delay interval, two consecutive time intervals without valid signals are selected on both sides of the positive delay signal window, respectively, as the first noise window and the second noise window corresponding to the positive delay; within the negative delay interval, two consecutive time intervals without valid signals are selected on both sides of the negative delay signal window, respectively, as the first noise window and the second noise window corresponding to the negative delay.

6. The seismic interferometric signal enhancement method based on an adaptive superposition strategy according to claim 1, characterized in that, The cross-spectral matrix refers to a two-dimensional matrix formed by arranging the cross-spectral results of all cross-correlation function pairs at frequency points according to the rule of pairing corresponding signals in rows and columns, using the filtered cross-correlation function data segments as indices.

7. The seismic interferometric signal enhancement method based on an adaptive superposition strategy according to claim 1, characterized in that, Filtered signal The calculation formula is: ， In the formula, The parameter used to adjust the roughness of the filter has a value greater than 0. This is a data segment representing the cross-correlation function in the frequency domain. To delay the transmission time, , For the index of stations, For time.

8. A seismic interferometric signal enhancement system based on an adaptive stacking strategy, employing the seismic interferometric signal enhancement method based on an adaptive stacking strategy as described in any one of claims 1-7, characterized in that, include: The cross-correlation calculation module is used to convert the seismic interferometric time-domain signal to the frequency domain and calculate the cross-correlation function data segment; The quality screening module is used to screen cross-correlation function data segments using the root mean square ratio screening superposition method; The cross-spectral analysis module is used to construct the cross-spectral matrix of the filtered data and design adaptive filters; The filtering module is used to filter the cross-correlation function data segment using an adaptive filter to obtain the filtered signal. The time-frequency analysis module is used to perform time-frequency decomposition on the filtered signal and extract instantaneous phase information in the frequency domain. The weight construction module is used to construct time-frequency domain coherent weights based on instantaneous phase information. The weighted superposition module is used to weight and superimpose the filtered signal using time-frequency domain coherent weights to obtain the superimposed time-frequency spectrum. The inverse transform module is used to inversely transform the superimposed time spectrum back to the time domain to obtain the final time-domain superposition result.

Images