A method and system for reconstructing and purifying seismic noise cross-correlation functions

Abstract

This invention discloses a method and system for reconstructing and purifying the noise cross-correlation function of seismic background noise, comprising the following steps: acquiring the original noise signal; utilizing the sparsity characteristics of the noise cross-correlation function of seismic background noise in the effective frequency band, and combining compressed sensing theory to perform non-uniform downsampling on the original noise signal to obtain undersampled data; employing a fast iterative shrinking threshold algorithm to reconstruct the full-band noise cross-correlation function from the undersampled data with high precision, obtaining the reconstructed noise cross-correlation function; performing frequency-wavenumber transformation on the reconstructed noise cross-correlation function, and using an FK filter mask to achieve three-component separation, obtaining the separated effective signal in the FK domain; performing an inverse transformation on the separated effective signal in the FK domain, and finally outputting the purified noise cross-correlation function.

Description

Technical Field

[0001] This invention belongs to the field of signal processing and geophysical inversion technology, specifically relating to a method and system for reconstructing and purifying seismic noise cross-correlation functions. Background Technology

[0002] Seismic background noise imaging technology is widely used in crustal structure inversion, resource exploration and other fields. Cross-correlation technology, represented by noise cross-correlation function, is an important part of the passive geophysical inversion field. Its calculation accuracy has become a key factor restricting the imaging accuracy, and thus affecting the correctness of the final geophysical inversion results.

[0003] Currently, traditional methods require uniform sampling across the entire frequency band to satisfy Nyquist's theorem, resulting in massive data volumes and high computational complexity. Methods for processing noise cross-correlation functions, such as linear Radon transform, all rely on dense sampling and have processing blind spots. High-frequency industrial noise signals, volume wave interference, and effective surface waves are aliased in the time-frequency domain, leading to effective signal loss using traditional filtering methods. These drawbacks are particularly pronounced when processing high-frequency signals from shallow ground. Furthermore, large-scale array data reconstruction is time-consuming and fails to meet real-time monitoring requirements, resulting in insufficient computational efficiency.

[0004] Existing methods for processing noise cross-correlation functions are rarely global, efficient, or geared towards signal inversion requirements. This often overlooked processing of the raw dataset represents a hidden technical challenge in seismic background noise imaging. Summary of the Invention

[0005] This invention aims to address the shortcomings of existing technologies and provides the following solutions:

[0006] A method for reconstructing and purifying the noise cross-correlation function of earthquake background noise, comprising the following steps:

[0007] The original noise signal is acquired, and the sparsity of the noise cross-correlation function of the earthquake background noise in the effective frequency band is utilized. Combined with compressed sensing theory, the original noise signal is non-uniformly downsampled to obtain undersampled data.

[0008] A fast iterative threshold shrinkage algorithm is used to reconstruct the full-band noise cross-correlation function from the undersampled data with high precision, resulting in the reconstructed noise cross-correlation function.

[0009] The reconstructed noise cross-correlation function is subjected to frequency-wavenumber transformation, and the three components are separated using an FK filter mask to obtain the separated effective signal in the FK domain.

[0010] The separated effective signal in the FK domain is inversely transformed to finally output the purified noise cross-correlation function.

[0011] Preferably, the data of the noise cross-correlation function includes: the latitude and longitude of the source station, the latitude and longitude of the receiving station, and the distance between the two stations.

[0012] Preferably, the method for performing the non-uniform downsampling includes:

[0013] Poisson disk random sampling or adaptive sampling based on signal energy distribution is adopted, and compressed sensing theory is introduced to perform non-uniform downsampling through the low coherence of the observation matrix and sparse basis to obtain undersampled data.

[0014] Preferably, the method for obtaining the reconstructed noise cross-correlation function includes:

[0015] Calculate the undersampled data:

[0016]

[0017] in, y This indicates undersampled data. x The time series representing the original noise cross-correlation function. Represents the observation matrix. n Indicates noise;

[0018] Based on the undersampled data, calculate the FISTA iterative optimization:

[0019]

[0020] in, F ( x () represents the objective function to be optimized. λ Represents the regularization parameter. Represents a sparse transformation basis;

[0021] Nesterov acceleration is performed during the calculation of the FISTA iterative optimization:

[0022]

[0023] in, x k Representation Algorithm Number k The sparse solution vector obtained after step iteration, prox λ|| .||1 This represents the soft threshold operator. z k Indicates the first k Auxiliary variables in the next iteration α Indicates the step size. t k This represents the sequence of time parameters controlling Nesterov momentum acceleration in the FISTA algorithm;

[0024] The reconstructed noise cross-correlation function is obtained by alternating iterative Nesterov accelerated gradient descent and soft threshold shrinkage.

[0025] Preferably, the effective signal in the separated FK domain includes: effective surface wave signal, high-frequency noise, and in-phase delayed volume wave at time 0.

[0026] The present invention also provides a system for reconstructing and purifying the noise cross-correlation function of earthquake background noise. The system applies the above-mentioned method and includes: a sampling module, a reconstruction module, a signal separation module, and an inverse transformation module.

[0027] The sampling module is used to acquire the original noise signal. It utilizes the sparsity of the noise cross-correlation function of the earthquake background noise in the effective frequency band and combines compressed sensing theory to perform non-uniform downsampling on the original noise signal to obtain undersampled data.

[0028] The reconstruction module is used to reconstruct the full-band noise cross-correlation function from the undersampled data with high precision using a fast iterative shrinking threshold algorithm, and obtain the reconstructed noise cross-correlation function.

[0029] The signal separation module is used to perform frequency-wavenumber transformation on the reconstructed noise cross-correlation function and use an FK filter mask to achieve three-component separation, thereby obtaining the separated effective signal in the FK domain.

[0030] The inverse transformation module is used to perform an inverse transformation on the separated effective signal in the FK domain, and finally outputs the purified noise cross-correlation function.

[0031] Preferably, the data of the noise cross-correlation function includes: the latitude and longitude of the source station, the latitude and longitude of the receiving station, and the distance between the two stations.

[0032] Preferably, the method for performing the non-uniform downsampling includes:

[0033] Poisson disk random sampling or adaptive sampling based on signal energy distribution is adopted, and compressed sensing theory is introduced to perform non-uniform downsampling through the low coherence of the observation matrix and sparse basis to obtain undersampled data.

[0034] Preferably, the method for obtaining the reconstructed noise cross-correlation function includes:

[0035] Calculate the undersampled data:

[0036]

[0037] in, y This indicates undersampled data. x The time series representing the original noise cross-correlation function. Represents the observation matrix.n Indicates noise;

[0038] Based on the undersampled data, calculate the FISTA iterative optimization:

[0039]

[0040] in, F ( x () represents the objective function to be optimized. λ Represents the regularization parameter. Represents a sparse transformation basis;

[0041] Nesterov acceleration is performed during the calculation of the FISTA iterative optimization:

[0042]

[0043] in, x k Representation Algorithm Number k The sparse solution vector obtained after step iteration, prox λ|| .||1 This represents the soft threshold operator. z k Indicates the first k Auxiliary variables in the next iteration α Indicates the step size. t k This represents the sequence of time parameters controlling Nesterov momentum acceleration in the FISTA algorithm;

[0044] The reconstructed noise cross-correlation function is obtained by alternating iterative Nesterov accelerated gradient descent and soft threshold shrinkage.

[0045] Preferably, the effective signal in the separated FK domain includes: effective surface wave signal, high-frequency noise, and in-phase delayed volume wave at time 0.

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

[0047] (1) The noise cross-correlation function of this invention has a high noise suppression capability. The fk joint domain notch filtering can specifically eliminate industrial frequency interference, surface wave pollution, and strong interference from volume wave effects, resulting in an extremely high signal-to-noise ratio. (2) This invention achieves a high signal-to-noise ratio at a 70% sampling rate through the improved FISTA algorithm, which is significantly better than traditional methods. The phase error is controlled within an acceptable range, and noise energy is significantly attenuated while retaining effective signal energy. (3) This invention provides efficient and high-precision reliable technical support for fields such as cross-correlation high-precision calculation, passive geophysical inversion such as seismic background noise imaging technology, and resource exploration. Attached Figure Description

[0048] To more clearly illustrate the technical solution of the present invention, the drawings used in the embodiments are briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0049] Figure 1 This is a schematic diagram of the method flow according to an embodiment of the present invention;

[0050] Figure 2 This is a comparison diagram of the reconstructed noise cross-correlation signal waveform, energy before and after reconstruction, and error distribution according to an embodiment of the present invention. Figure 2 (a) in the diagram represents the cross-correlation function of the original input noise. Figure 2 (b) in the diagram represents the signal region. Figure 2 (c) in the diagram represents the energy difference map in the FK domain. Figure 2 In the diagram, (d) represents the reconstructed noise cross-correlation function. Figure 2 In the diagram, (e) represents the energy difference map before and after reconstruction. Figure 2 In the diagram, (f) represents the error distribution plot;

[0051] Figure 3 The following are examples of the FK domain spectra before and after reconstruction of the noise cross-correlation function in this embodiment of the invention, the components retained in the reconstructed FK domain, and the filter mask diagram. Figure 3 In the diagram, (a) represents the original noise cross-correlation function FK spectrum. Figure 3 (b) in the image represents the FK spectrum after CS reconstruction. Figure 3 (c) in the image represents the spectrum after FK denoising. Figure 3 In this context, (d) represents the FK filter mask;

[0052] Figure 4 This diagram illustrates the power spectral density comparison and frequency suppression effect detection of the original noise cross-correlation function, the noise cross-correlation function reconstructed by compressed sensing, and the noise cross-correlation function after FK filtering at inter-station spacings of 10.02 km, 20.54 km, and 31.05 km, according to an embodiment of the present invention. Figure 4 (a) shows a schematic diagram comparing the power spectral density of the three noise cross-correlation functions at a distance of 10.02 km between the stations. Figure 4 (b) shows a schematic diagram illustrating the frequency suppression effect of the three noise cross-correlation functions at a station spacing of 10.02 km. Figure 4 (c) represents a schematic diagram comparing the power spectral density of the three noise cross-correlation functions at a distance of 20.54 km between the stations. Figure 4(d) represents a schematic diagram illustrating the frequency suppression effect of the three noise cross-correlation functions at a station spacing of 20.54 km. Figure 4 (e) represents a schematic diagram comparing the power spectral density of the three noise cross-correlation functions at a distance of 31.05 km between the stations. Figure 4 In the diagram (f), the frequency suppression effect of the three noise cross-correlation functions at a distance of 31.05 km between stations is shown.

[0053] Figure 5 This diagram shows a comparison of the waveform amplitudes of the original noise cross-correlation function, the noise cross-correlation function reconstructed by compressed sensing, and the noise cross-correlation function after FK filtering at inter-station spacings of 2.5 km, 4.01 km, and 12.02 km, along with magnified detail window diagrams. Figure 5 (a) shows a schematic diagram comparing the waveform amplitudes of the three noise cross-correlation functions at a distance of 2.5 km between the stations. Figure 5 (b) shows a magnified view of the detail window of the three noise cross-correlation functions at a distance of 2.5 km between stations. Figure 5 (c) shows a schematic diagram comparing the waveform amplitudes of the three noise cross-correlation functions at a distance of 4.01 km between the stations. Figure 5 (d) represents a magnified view of the detail window of the three noise cross-correlation functions at a distance of 4.01 km between stations. Figure 5 (e) represents a schematic diagram comparing the waveform amplitudes of the three noise cross-correlation functions at a distance of 12.02 km between the stations. Figure 5 In the diagram (f), we see a magnified view of the detail window of the three noise cross-correlation functions at a spacing of 12.02 km between stations.

[0054] Figure 6 This diagram illustrates a comparison of the waveform amplitudes and power spectral densities of the original noise cross-correlation function, the noise cross-correlation function reconstructed by compressed sensing, and the noise cross-correlation function after FK filtering at inter-station spacings of 2.5 km, 4.01 km, and 12.02 km, according to an embodiment of the present invention. Figure 6 (a) shows a schematic diagram comparing the waveform amplitudes of the three noise cross-correlation functions at a distance of 2.5 km between the stations. Figure 6 (b) shows a schematic diagram comparing the power spectral density of the three noise cross-correlation functions at a distance of 2.5 km between the stations. Figure 6 (c) shows a schematic diagram comparing the waveform amplitudes of the three noise cross-correlation functions at a distance of 4.01 km between the stations. Figure 6 (d) represents a schematic diagram comparing the power spectral density of the three noise cross-correlation functions at a distance of 4.01 km between the stations. Figure 6 (e) represents a schematic diagram comparing the waveform amplitudes of the three noise cross-correlation functions at a distance of 12.02 km between the stations. Figure 6 In the diagram (f), we see a comparison of the power spectral density of the three noise cross-correlation functions at a distance of 12.02 km between the stations. Detailed Implementation

[0055] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0056] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0057] Example 1

[0058] In this embodiment, as Figure 1 As shown, a method for reconstructing and purifying the noise cross-correlation function of earthquake background noise includes the following steps:

[0059] The original noise signal is collected, and the sparsity of the noise cross-correlation function of the earthquake background noise in the effective frequency band is utilized. Combined with compressed sensing theory, the original noise signal is non-uniformly downsampled to obtain undersampled data.

[0060] The noise cross-correlation function data includes: the latitude and longitude of the source station, the latitude and longitude of the receiving station, and the distance between the two stations. Methods for non-uniform downsampling include: using Poisson disk random sampling or adaptive sampling based on signal energy distribution, while simultaneously introducing compressed sensing theory to perform non-uniform downsampling through the low coherence of the observation matrix and sparse basis, resulting in undersampled data.

[0061] In this embodiment, the sparse characteristics of the noise cross-correlation function in the effective frequency band are utilized, combined with compressed sensing theory, to perform non-uniform downsampling of the original signal. Within a specific frequency band, the effective signal energy is concentrated in a few frequency components, while the energy in other frequency bands approaches zero. Its mathematical representation lies in the fact that, after Fourier transform or wavelet transform, the coefficient distribution of the signal in the frequency domain exhibits a sparse pattern of "a few large values ​​+ a large number of small values," which conforms to the definition of sparse signals in compressed sensing theory. L(0-norm minimization). In the FK domain noise cross-correlation function, the effective surface wave signal exhibits low wavenumber and low-frequency energy clusters, which are naturally separable from high-frequency noise (high wavenumber diffuse distribution) and volume waves (high wavenumber linear events) in the space-frequency joint domain, providing a physical basis for filtering. The sparsity assumption is that the signal has only a few non-zero coefficients in a certain transform domain (such as the frequency domain, wavelet domain, or FK domain), with the remaining coefficients close to zero. The sparsity characterization is performed by performing an FK transform or wavelet transform on the noise cross-correlation function to confirm the sparse basis.

[0062] Sparsity allows the noise cross-correlation function to retain key information in the effective frequency band through non-uniform sampling, overcoming the Nyquist sampling rate limit. The time series of the noise cross-correlation function is sampled using a Poisson disk or random intervals to avoid spectral aliasing caused by periodic repetition.

[0063] A fast iterative threshold shrinkage algorithm is used to reconstruct the full-band noise cross-correlation function from undersampled data with high accuracy, resulting in the reconstructed noise cross-correlation function.

[0064] In this embodiment, based on the sparsity of the noise cross-correlation function in the effective frequency band, the original signal is non-uniformly downsampled (e.g., retaining 70% of the data points); then, through FISTA iterative optimization, combined with gradient descent and soft thresholding shrinkage operations, the high-resolution signal is gradually approximated from the undersampled data—gradient descent ensures data fitting, and soft thresholding forces sparsity to suppress noise.

[0065] Specifically, methods for obtaining the reconstructed noise cross-correlation function include: calculating the undersampled data:

[0066]

[0067] in, y This indicates undersampled data. x The time series representing the original noise cross-correlation function. Represents the observation matrix. n Noise is represented; based on undersampled data, FISTA iterative optimization is calculated. Specifically, FISTA is an iterative algorithm for solving sparse optimization problems, used to solve regularization problems of the following form:

[0068]

[0069] in, F ( x Let ) represent the objective function to be optimized, which is expressed by the data fitting term and the sparse regularization term. λ This represents the regularization parameter, which controls the sparsity intensity. Denotes a sparse transformation basis. for f ( x The quadratic term of ) for g ( x ) sparse terms, f ( x ) represents a differentiable term. g ( x ) represents a nondifferentiable term; the number of iterations in FISTA refers to the number of loops in which the algorithm updates the sparse solution vector. Each iteration includes: (1) Gradient descent step: gradient update, calculating the gradient direction of the current solution. , T Represents the transpose of the matrix; updates the solution along the negative gradient direction of the loss function. Where L is the Lipschitz constant, taken as... The spectral norm can be estimated in practice using the power iteration method. (2) Shrinking threshold step: Apply L 1 constraint ,in, The threshold function is used. Sparsity is achieved by applying a soft threshold operator. (3) Momentum acceleration step: Momentum acceleration, introducing auxiliary variables. y k Achieving supergradient descent During the FISTA iterative optimization process, Nesterov acceleration is performed, and the step size is adjusted using Nesterov acceleration techniques. Introducing a momentum term into the gradient descent iterations improves the convergence speed from... Upgraded to , k Let be the number of iterations, expressed as:

[0070] in, x k Representation Algorithm Number k The sparse solution vector obtained after one iteration ensures the sparsity of the solution through a soft thresholding operation (proximal mapping); prox λ|| .||1 This represents the soft thresholding operator, which shrinks each element of the input vector toward zero; z k Indicates the first k The auxiliary variables (or extrapolation points) in the next iteration, and the main variables x k Together they form an iterative sequence; α Indicates the step size, usually taken as This can be optimized through line search; t k This represents the time parameter sequence controlling Nesterov momentum acceleration in the FISTA algorithm, with initial values ​​typically set to... t 0=1; Under high signal-to-noise ratio, The convergence threshold is The solution is obtained by alternating iterative steps of Nesterov accelerated gradient descent and soft threshold shrinkage. After 200 iterations, a sparse solution x200 is finally output. When the iteration reaches 200, the following condition is typically met: At this point, increasing the number of iterations effectively improves the quality of the solution. 200 iterations are sufficient to stably converge to Gaussian white noise. However, impulse noise requires 300+ iterations to overcome local minima. FISTA's... The convergence rate makes 200 iterations equivalent to 800-1000 iterations of ISTA (Iterative Shrinkage-Thresholding Algorithm). Finally, the reconstructed noise cross-correlation function is obtained.

[0071] The reconstructed noise cross-correlation function is subjected to frequency-wavenumber transformation, and the three components are separated using an FK filter mask to obtain the separated effective signal in the FK domain.

[0072] The effective signal in the separated FK domain includes: effective surface wave signal, high-frequency noise, and in-phase delayed volume wave at time 0.

[0073] In this embodiment, the cross-correlation function of time-domain noise is transformed to the FK domain. Utilizing the distribution differences of surface waves (low frequency, low wavenumber), high-frequency noise (wideband, high wavenumber), and volume waves (in-phase axis at time 0, exhibiting a specific linear trajectory) in the frequency-wavenumber space, an adaptive filtering mask is constructed. This mask retains effective surface waves through a low-wavenumber passband and suppresses high-frequency noise through a high-wavenumber bandstop. A notch filter is designed to address the wavenumber characteristics of volume wave interference. Finally, the separated time-domain signal is obtained through inverse transformation, achieving signal-to-noise ratio improvement and volume wave interference elimination, providing a clean data foundation for subsequent wave velocity measurements.

[0074] The effective signal in the separated FK domain is inversely transformed to finally output the purified noise cross-correlation function.

[0075] In this embodiment, the surface wave dominant signal (preserving low-frequency and low-wavenumber components) after FK filtering mask processing is subjected to inverse frequency-wavenumber transformation to convert it from the wavenumber-frequency joint domain back to the time-space domain. During this process, special attention should be paid to the accurate reconstruction of phase information, and windowed Fourier inverse transform or fast FK inversion algorithm is used to ensure waveform integrity. The final output purified noise cross-correlation function has three main characteristics: (1) the effective surface wave dispersion characteristics are completely preserved, (2) the high-frequency noise components are significantly attenuated, and (3) the in-phase axial body wave interference at time 0 is basically eliminated, providing high-quality input data that meets the requirements of Green's function theory for subsequent crustal structure inversion.

[0076] Example 2

[0077] In this embodiment, a noise cross-correlation function reconstruction and purification system for earthquake background noise includes: a sampling module, a reconstruction module, a signal separation module, and an inverse transformation module.

[0078] The sampling module is used to acquire the original noise signal. It utilizes the sparsity of the noise cross-correlation function of the earthquake background noise in the effective frequency band and combines compressed sensing theory to perform non-uniform downsampling on the original noise signal to obtain undersampled data.

[0079] The noise cross-correlation function data includes: the latitude and longitude of the source station, the latitude and longitude of the receiving station, and the distance between the two stations.

[0080] Methods for non-uniform downsampling include: using Poisson disk random sampling or adaptive sampling based on signal energy distribution, while introducing compressed sensing theory, and using the low coherence of the observation matrix and sparse basis to perform non-uniform downsampling to obtain undersampled data.

[0081] The reconstruction module is used to reconstruct the full-band noise cross-correlation function from undersampled data with high precision using a fast iterative shrinking threshold algorithm, and obtain the reconstructed noise cross-correlation function.

[0082] Methods for obtaining the reconstructed noise cross-correlation function include: calculating undersampled data:

[0083]

[0084] in, y This indicates undersampled data. x The time series representing the original noise cross-correlation function. Represents the observation matrix. n Representing noise; calculating FISTA iterative optimization based on undersampled data:

[0085]

[0086] in, F ( x () represents the objective function to be optimized. λ Represents the regularization parameter. Represents the sparse transformation basis; Nesterov acceleration is performed during the computation of FISTA iterative optimization:

[0087]

[0088] in, x k Representation Algorithm Number k The sparse solution vector obtained after step iteration, prox λ|| .||1 This represents the soft threshold operator. z k Indicates the firstk Auxiliary variables in the next iteration α Indicates the step size. t k This represents the time parameter sequence controlling Nesterov momentum acceleration in the FISTA algorithm; the reconstructed noise cross-correlation function is obtained by iteratively solving Nesterov accelerated gradient descent and soft threshold contraction.

[0089] The signal separation module is used to perform frequency-wavenumber transformation on the reconstructed noise cross-correlation function and to achieve three-component separation using an FK filter mask to obtain the separated effective signal in the FK domain.

[0090] The effective signal in the separated FK domain includes: effective surface wave signal, high-frequency noise, and in-phase delayed volume wave at time 0.

[0091] The inverse transform module is used to perform an inverse transform on the separated effective signal in the FK domain, and finally outputs the purified noise cross-correlation function.

[0092] Example 3

[0093] In this embodiment, a device for reconstructing and purifying noise cross-correlation functions based on compressed sensing fast iterative shrinkage threshold iteration and frequency-wavenumber division is also provided, including:

[0094] This toolkit is designed for reading, writing, and standardizing seismic SAC format data. It handles batch loading, distance sorting, format conversion, and result output of seismic data. Standardized data loading integrates scattered SAC files into a unified time-distance panel format, facilitating subsequent FK domain analysis or other seismic signal processing. Spatial relationship processing automatically calculates epicentral distances and sorts them by distance, establishing a spatially ordered data structure. All waveform data is uniformly truncated to the same length to ensure dimensionality matching in subsequent processing and guarantee data consistency.

[0095] The seismic signal processing device provides a general toolchain module covering five major functions: log management, performance monitoring, data preprocessing, signal analysis and quality assessment, supporting the complete process of seismic data denoising, reconstruction and feature extraction.

[0096] A seismic data denoising device based on the FK domain is used for noise suppression and effective wave enhancement in seismic signal processing. Seismic data denoising eliminates linear noise (such as surface waves and multiples) and random noise in seismic records through FK domain filtering, preserving effective reflection signals. Non-uniform data normalization includes processing non-uniformly distributed data acquired in the field, achieving uniform spatial intervals through FK analysis via resampling. Velocity-selective filtering involves designing velocity bandpass filters based on differences in seismic wave propagation velocities.

[0097] Non-uniform sampling devices reduce data volume or focus on key signal features through probabilistic sampling. Data compression reduces data volume through proportional downsampling, suitable for optimizing seismic data storage or transmission. Energy-oriented processing prioritizes the retention of high-energy segments and suppresses low-energy noise. The algorithm provides a baseline for energy-sensing sampling, using uniform sampling as a reference.

[0098] A compressed sensing reconstruction framework based on a fast iterative shrinking threshold algorithm is proposed for sparse recovery of one-dimensional signals and extended to the processing of two-dimensional seismic gathers. It recovers complete signals from non-uniformly sampled data, reduces data acquisition costs, and achieves efficient reconstruction by leveraging the sparsity assumption of Discrete Cosine Transform (DCT) and the fast convergence characteristics of FISTA.

[0099] Compressed sensing reconstruction and FK domain denoising device for recovering high-quality signals and suppressing noise from sparsely sampled seismic data. It recovers complete signals from non-uniformly sampled seismic data using compressed sensing technology, identifies and filters noise in the FK domain, retains effective signals, and integrates a one-stop solution for data loading, sampling, reconstruction, denoising, visualization and quality assessment.

[0100] The earthquake data visualization toolset provides comprehensive visualization support for the earthquake data processing workflow. It compares the spatiotemporal-frequency characteristics of the original data, compressed sampling reconstruction results, and denoised results through multimodal charts, assisting in algorithm performance evaluation and parameter optimization.

[0101] Example 4

[0102] In this embodiment, a method for reconstructing and refining noise cross-correlation functions based on a compressed sensing fast iterative threshold shrinkage algorithm and frequency-wavenumber division includes:

[0103] Step S1: Obtain the noise cross-correlation function and perform preprocessing to obtain noise cross-correlation function data containing surface wave signals;

[0104] Step S2: Integrate the scattered SAC files into a time-distance panel of a unified format to facilitate subsequent FK domain analysis or other seismic signal processing. Automatically calculate the epicentral distance and sort it by distance to establish a spatially ordered data structure. Truncate all waveform data to the same length to ensure dimensionality matching in subsequent processing. Extract a consistent time window and apply a cosine cone window to suppress edge effects.

[0105] Step S3: Employ square-dimensional units to enhance the differences in high-energy areas, converting energy distribution into probability density to ensure the rationality of irregular sampling. Reduce data volume through proportional downsampling, suitable for optimizing seismic data storage or transmission. Prioritize the retention of high-energy segments, suppress low-energy noise, and provide uniform sampling as a reference baseline for energy sensing sampling. Perform completely random uniform sampling based on the probability distribution of waveform energy values ​​(higher energy values ​​have a higher probability of being sampled).

[0106] Step S4: Perform FISTA-DCT reconstruction, where the objective function is:

[0107]

[0108] in, For the signal to be reconstructed, This represents the observation signal with missing data (unsampled points are set to 0). For diagonal mask matrix ( when i (Points are sampled) The DCT transformation matrix is... λ For sparse regularization parameters.

[0109] Based on input observed signal and mask M Initialize variables:

[0110]

[0111]

[0112]

[0113] in, This represents the Nesterov momentum auxiliary variable. It is the first k The extrapolation point (auxiliary variable) in the next iteration is used to accelerate convergence. It is a time-series parameter that controls the momentum weight, and is dynamically adjusted. The proportion of historical information being updated, balancing convergence speed and stability, and initial values. The step size (Lipschitz constant) is (The matrix spectral norm needs to be calculated or estimated.)

[0114] Based on the FISTA main iteration, for Perform gradient calculation:

[0115]

[0116] in, ;

[0117] Accelerated updates based on Nesterov:

[0118]

[0119] in, Indicates a temporary variable. It indicates that it is the first k Extrapolation point in the next iteration;

[0120] Proximal mapping: Element-wise soft thresholding function: ;

[0121] Update based on inverse transform:

[0122] ;

[0123] Momentum update:

[0124]

[0125] ;

[0126] The iteration is complete if the maximum number of iterations or the relative error condition is met.

[0127] .

[0128] Step S5: Automatically detect sparse arrays and switch to conservative filtering mode. Support standardized processing of non-uniform spatial sampling data. Add a two-dimensional cosine cone window to reduce spectral leakage, pad with zeros to powers of 2, calculate the spectrum using FFT2, center it, multiply the velocity sector mask with the spectrum, restore the time-domain signal using inverse FFT, interpolate the uniform gather back to the original non-uniform spatial location, and record key parameters such as spatial sampling rate, Nyquist wavenumber, and velocity resolution.

[0129] Step S6: Display the time-distance profiles of the original / sampled / reconstructed / denoised data side by side, compare the frequency-wavenumber domain energy distribution of the original / reconstructed / denoised data, extract single-channel waveforms according to the distance ratio, compare the three-stage signal morphology, extract the amplitude of all channels at fixed time points, and analyze spatial consistency.

[0130] In reconstructing data, Figure 2 The reconstruction noise cross-correlation signal comparison, energy difference before and after, and error distribution diagram (RawPanel is the original input noise cross-correlation function, FK Denoised Panel is the reconstructed new noise cross-correlation function. The reconstructed new noise cross-correlation function of FK Denoised Panel effectively eliminates the volume wave effect at time 0 in the noise cross-correlation function, and weakens a large amount of high-frequency noise while retaining the effective surface wave signal).

[0131] The main noise component was eliminated in the FK domain using a filtering mask. Figure 3 ), Figure 4 The power spectral density (vertical axis represents power spectral density) of the waveforms at inter-station spacings of 10.02 km, 20.54 km, and 31.05 km is compared for the original noise cross-correlation function, the noise cross-correlation function reconstructed by compressed sensing, and the noise cross-correlation function after FK filtering, and the frequency suppression effect is detected. After reconstruction, the power of the high-frequency part is reduced while retaining the effective signal frequency.

[0132] Figure 5 The waveform amplitudes of the original noise cross-correlation function, the noise cross-correlation function reconstructed by compressed sensing, and the noise cross-correlation function after FK filtering are compared at station spacings of 2.5km, 4.01km, and 12.02km, and detail window magnification is also provided. Figure 6 For comparative analysis of the time-domain waveforms and power spectra before and after reconstruction: Waveform comparison at a station spacing of 2.5km (blue solid line: original noise cross-correlation function; green dashed line: reconstructed waveform after FK domain filtering; orange solid line: CS reconstructed waveform); local time window magnification of the waveform at 2.5km; comparison of the original noise cross-correlation function, the power spectrum of the compressed sensing reconstructed waveform, and the FK domain filtered reconstructed waveform at 2.5km; Waveform comparison at a station spacing of 4.01km (blue solid line: original noise cross-correlation function; green dashed line: reconstructed waveform after FK domain filtering; orange solid line: compressed sensing reconstructed waveform); local time window magnification of the waveform at 4.01km; comparison of the original noise cross-correlation function, the compressed sensing reconstructed waveform, and the FK domain filtered reconstructed waveform at 4.01km; Local time window magnification of the waveform at 12.02km; comparison of the original noise cross-correlation function, the compressed sensing reconstructed waveform, and the FK domain filtered reconstructed waveform at 12.02km.

[0133] Figure 3 , Figure 6 The processing effectiveness of compressed sensing reconstruction and subsequent FK domain filtering was systematically evaluated across three key dimensions: FK spectrum evolution, time-domain waveform optimization, and spectral energy redistribution. In the FK domain analysis, the selective filtering mask effectively eliminated dominant noise energy while preserving coherent signal components. (At short station spacing...) Figure 6 At a distance of 2.5 km, coherent noise near zero delay is significantly suppressed without altering the frequency components and amplitude characteristics of the main signal; at larger station spacing ( Figure 6 (12.02km), the filtering process efficiently attenuates high-frequency noise energy while completely preserving the time-frequency characteristics of the main surface wave energy. Figure 6(4.01 km). Power spectral density analysis further shows that the reconstruction process preserves low-frequency signal components (<3 Hz) with minimal loss, while significantly suppressing the high-frequency band (>3 Hz), which is dominated by incoherent noise. These results collectively confirm that the reconstruction based on compressed sensing and FK domain filtering achieves effective noise suppression while ensuring spectral fidelity and signal coherence. This provides a more physically consistent dataset for subsequent tomographic inversion, thereby ensuring higher stability and resolution of the inversion speed model.

[0134] Example 5

[0135] The program developed based on this invention provides complete compressed sensing fast iterative shrinkage threshold iteration and frequency-wavenumber division noise cross-correlation function reconstruction functions, specifically:

[0136] The data processing flow includes input and output, with key steps including data loading: reading and sorting the SAC file to generate a data panel of uniform length. Energy sensing sampling: compressing and sampling the raw signal proportionally.

[0137] Compressed Sensing Reconstruction: The FISTA algorithm is iterated 200 times to reconstruct missing data. FK Denoising: Noise is filtered by setting a velocity band, preserving the valid signal. Visualization and Reporting: Analytical charts such as spatiotemporal domain, FK spectrum, and profile plots are generated and saved as PDF reports. Core algorithms and parameters include: Compressed Sensing Regularization parameter: lambda_reg=0.01, Iteration count: iters=200

[0138] Reconstruction results. FK domain denoising: velocity resolution. Wavenumber domain mask: generated via velocity bands, covering 11.5% of the passband.

[0139] Energy after denoising: The energy of the ROI region decreased from the original 4879.575 to 901.558. Full-process visualization: Provides multi-dimensional spatiotemporal and frequency domain analysis tools.

[0140] Input parameters include: --erase_roi, specifying the time-distance rectangular region [t0,t1,d0,d1] (seconds, kilometers), used to suppress signals in specific areas (such as direct wave interference). --edge_t, the soft-edge transition width of the time axis (seconds), to avoid the Gibbs effect of rectangular filtering. --edge_d, the soft-edge transition width of the distance axis (kilometers). --in (input_dir), the input directory path (must contain SAC format seismic data). --out (output_dir), the output directory path (the location where the results files are stored). --ratio, the compressed sensing sampling ratio (0-1), controlling the data downsampling rate. --lambda (lambda_reg), the compressed sensing regularization parameter, balancing data fidelity and sparsity (the larger the value, the smoother the signal). --iters, the maximum number of iterations for the FISTA algorithm. --bands, the FK filter velocity bands (km / s), in nested list format, used to retain signals within a specified velocity range. --f_notch, the frequency notch threshold (Hz), to remove low-frequency noise (e.g., <0.1Hz). --k_notch, wavenumber notch threshold (1 / km), suppresses near-field interference. --smooth, FK mask Gaussian smoothing coefficient (pixels), avoids artifacts caused by sharp cutoff. --seed, random seed, ensures the sampling process is reproducible. --profile_locs, waveform comparison profile positions (relative distance ratio), used for visualization analysis.

[0141] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Various modifications and improvements made to the technical solutions of the present invention by those skilled in the art without departing from the spirit of the present invention should fall within the protection scope defined by the claims of the present invention.

Claims

1. A method of noise cross-correlation function reconstruction and purification of seismic ambient noise, characterized in that, Includes the following steps: The original noise signal is acquired, and the sparsity of the noise cross-correlation function of the earthquake background noise in the effective frequency band is utilized. Combined with compressed sensing theory, the original noise signal is non-uniformly downsampled to obtain undersampled data. A fast iterative threshold shrinkage algorithm is used to reconstruct the full-band noise cross-correlation function from the undersampled data with high precision, resulting in the reconstructed noise cross-correlation function. The reconstructed noise cross-correlation function is subjected to frequency-wavenumber transformation, and the three components are separated using an FK filter mask to obtain the separated effective signal in the FK domain. The separated effective signal in the FK domain is inversely transformed to finally output the purified noise cross-correlation function.

2. The method of claim 1, wherein, The noise cross-correlation function data includes: the latitude and longitude of the source station, the latitude and longitude of the receiving station, and the distance between the two stations.

3. The method of claim 1, wherein, The method for performing the non-uniform downsampling includes: Poisson disk random sampling or adaptive sampling based on signal energy distribution is adopted, and compressed sensing theory is introduced to perform non-uniform downsampling through the low coherence of the observation matrix and sparse basis to obtain undersampled data.

4. The method of claim 1, wherein, The method for obtaining the reconstructed noise cross-correlation function includes: Calculate the undersampled data: wherein, y denotes the under-sampled data, x denotes the time series of the original noise cross-correlation function, denotes the observation matrix, n denotes the noise; Based on the undersampled data, calculate the FISTA iterative optimization: in, F ( x () represents the objective function to be optimized. λ Represents the regularization parameter. Denotes a sparse transformation basis; Nesterov acceleration is performed during the calculation of the FISTA iterative optimization: in, x k Representation Algorithm Number k The sparse solution vector obtained after step iteration, prox λ|| .||1 This represents the soft threshold operator. z k Indicates the first k Auxiliary variables in the next iteration α Indicates the step size. t k This represents the sequence of time parameters controlling Nesterov momentum acceleration in the FISTA algorithm; The reconstructed noise cross-correlation function is obtained by alternating iterative Nesterov accelerated gradient descent and soft threshold shrinkage.

5. The method for reconstructing and purifying the noise cross-correlation function of earthquake background noise according to claim 1, characterized in that, The effective signal in the separated FK domain includes: effective surface wave signal, high-frequency noise, and in-phase delayed volume wave at time 0.

6. A system for reconstructing and purifying the noise cross-correlation function of earthquake background noise, wherein the system applies the method described in any one of claims 1-5, characterized in that, include: Sampling module, reconstruction module, signal separation module, and inverse transform module; The sampling module is used to acquire the original noise signal. It utilizes the sparsity of the noise cross-correlation function of the earthquake background noise in the effective frequency band and combines compressed sensing theory to perform non-uniform downsampling on the original noise signal to obtain undersampled data. The reconstruction module is used to reconstruct the full-band noise cross-correlation function from the undersampled data with high precision using a fast iterative shrinking threshold algorithm, and obtain the reconstructed noise cross-correlation function. The signal separation module is used to perform frequency-wavenumber transformation on the reconstructed noise cross-correlation function and use an FK filter mask to achieve three-component separation, thereby obtaining the separated effective signal in the FK domain. The inverse transformation module is used to perform an inverse transformation on the separated effective signal in the FK domain, and finally outputs the purified noise cross-correlation function.

7. The system for reconstructing and purifying the noise cross-correlation function of earthquake background noise according to claim 6, characterized in that, The noise cross-correlation function data includes: the latitude and longitude of the source station, the latitude and longitude of the receiving station, and the distance between the two stations.

8. The system for reconstructing and purifying the noise cross-correlation function of seismic background noise according to claim 6, characterized in that, The method for performing the non-uniform downsampling includes: Poisson disk random sampling or adaptive sampling based on signal energy distribution is adopted, and compressed sensing theory is introduced to perform non-uniform downsampling through the low coherence of the observation matrix and sparse basis to obtain undersampled data.

9. The system for reconstructing and purifying the noise cross-correlation function of earthquake background noise according to claim 6, characterized in that, The method for obtaining the reconstructed noise cross-correlation function includes: Calculate the undersampled data: in, y This indicates undersampled data. x The time series representing the original noise cross-correlation function. Represents the observation matrix. n Indicates noise; Based on the undersampled data, calculate the FISTA iterative optimization: in, F ( x () represents the objective function to be optimized. λ Represents the regularization parameter. Denotes a sparse transformation basis; Nesterov acceleration is performed during the calculation of the FISTA iterative optimization: in, x k Representation Algorithm Number k The sparse solution vector obtained after step iteration, prox λ|| .||1 This represents the soft threshold operator. z k Indicates the first k Auxiliary variables in the next iteration α Indicates the step size. t k This represents the sequence of time parameters controlling Nesterov momentum acceleration in the FISTA algorithm; The reconstructed noise cross-correlation function is obtained by alternating iterative Nesterov accelerated gradient descent and soft threshold shrinkage.

10. The system for reconstructing and purifying the noise cross-correlation function of seismic background noise according to claim 6, characterized in that, The effective signal in the separated FK domain includes: effective surface wave signal, high-frequency noise, and in-phase delayed volume wave at time 0.

Images