Seismic noise cross-correlation function reconstruction and purification method and system

Abstract

The invention discloses a method and system for reconstructing and purifying a noise cross-correlation function of seismic background noise, and the method comprises the following steps: collecting an original noise signal, carrying out the non-uniform downsampling of the original noise signal through employing the sparse characteristic of the noise cross-correlation function of the seismic background noise in an effective frequency band, and combining with a compressed sensing theory, and obtaining a non-uniform downsampling signal; obtaining under-sampling data; reconstructing a full-band noise cross-correlation function from the undersampled data with high precision by adopting a fast iterative shrinkage threshold algorithm to obtain a reconstructed noise cross-correlation function; performing frequency-wave number transformation on the reconstructed noise cross-correlation function, and realizing three-component separation by using an FK filtering mask to obtain a separated FK domain effective signal; and carrying out inverse transformation on the separated FK domain effective signal, and finally outputting a purified noise cross-correlation function.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of signal processing and geophysical inversion, and particularly relates to a method and system for reconstructing and purifying seismic noise cross-correlation functions. BACKGROUND

[0002] The cross-correlation technology represented by the noise cross-correlation function is an important part in the field of passive geophysical inversion, and its calculation accuracy is a key factor restricting the imaging accuracy, which further affects the correctness of the final result of subsequent geophysical inversion.

[0003] At present, the traditional method needs to uniformly sample all frequency bands to meet the Nyquist theorem, resulting in a large amount of data and high computational complexity. The noise cross-correlation function processing methods such as linear Radon transform all rely on dense sampling, and there are processing blind areas; the industrial noise high-frequency signal, body wave interference and effective surface wave are mixed in the time-frequency domain, and the traditional filtering method leads to the loss of effective signal, and these shortcomings are particularly prominent in processing shallow high-frequency signals; the reconstruction of large-scale array data takes a long time and is difficult to meet the real-time monitoring demand, and the computational efficiency is insufficient.

[0004] In the existing noise cross-correlation function processing, there are few processing methods that are global, efficient and oriented to signal inversion requirements. This easily overlooked raw data set processing is a hidden technical challenge in seismic ambient noise imaging. SUMMARY

[0005] The present application aims to solve the problems of the prior art, and provides the following scheme: A method for reconstructing and purifying noise cross-correlation functions of seismic ambient noise, comprising the following steps: Collecting original noise signals, using the sparse characteristics of noise cross-correlation functions of seismic ambient noise in the effective frequency band, and combining the compressed sensing theory to non-uniformly downsample the original noise signals to obtain undersampled data; Using a fast iterative shrinkage threshold algorithm to reconstruct the full-band noise cross-correlation function from the undersampled data with high accuracy to obtain the reconstructed noise cross-correlation function; Performing frequency-wavenumber transformation on the reconstructed noise cross-correlation function, and using FK filter mask to realize three-component separation to obtain separated FK domain effective signals; Performing inverse transformation on the separated FK domain effective signals to finally output the purified noise cross-correlation function.

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

[0007] Preferably, the method for performing the non-uniform decimation comprises: By using Poisson disc random sampling or adaptive sampling based on signal energy distribution, and by introducing the theory of compressive sensing, non-uniform decimation is performed through the low coherence of the observation matrix and the sparse basis to obtain the under-sampling data.

[0008] Preferably, the method for obtaining the reconstructed noise cross-correlation function comprises: The under-sampling data is calculated as: wherein, y represents the under-sampling data, x represents the time series of the original noise cross-correlation function, represents the observation matrix, n represents the noise; Based on the under-sampling data, the FISTA iterative optimization is calculated as: wherein, F represents the objective function to be optimized, x represents the regularization parameter, λ represents the sparse transform basis; During the calculation of the FISTA iterative optimization, Nesterov acceleration is performed as: wherein, k represents the sparse solution vector obtained after the algorithm iteration, x k represents the soft threshold operator, k represents the auxiliary variable in the iteration, z represents the step size, k k represents the time parameter sequence for controlling the Nesterov momentum acceleration in the FISTA algorithm; α The reconstructed noise cross-correlation function is obtained by solving through the alternating iteration of the Nesterov accelerated gradient descent and the soft threshold shrinkage. t Preferably, the effective signal in the FK domain after separation comprises: effective surface wave signal, high frequency noise and 0-time in-phase time delay body wave.

[0009] The application further provides a noise cross-correlation function reconstruction and purification system for seismic background noise, which applies the above method and comprises a sampling module, a reconstruction module, a signal separation module and an inverse transform module.

[0010] The application further provides a noise cross-correlation function reconstruction and purification system for seismic background noise, which applies the above method and comprises a sampling module, a reconstruction module, a signal separation module and an 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 to obtain 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.

[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: 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.

[0013] Preferably, 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. Represents 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 one iteration. denotes a soft threshold operator, z k denotes the first k auxiliary variable in the α denotes a step size, t k denotes a sequence of time parameters controlling Nesterov momentum acceleration in the FISTA algorithm; solved by alternating iterations of Nesterov accelerated gradient descent and soft thresholding shrinkage, to obtain the reconstructed noise cross-correlation function.

[0014] Preferably, the separated FK domain effective signal comprises: effective surface wave signal, high frequency noise and 0 time in-phase time delay body wave.

[0015] Compared with the prior art, the present application has the following beneficial effects: (1) The noise cross-correlation function of the present application has high noise suppression capability. The f-k joint domain wave filter can specifically eliminate industrial frequency point interference, surface wave pollution, and strong body wave effect interference, and has extremely strong signal-to-noise ratio. (2) The present application realizes high signal-to-noise ratio at 70% sampling rate by using the improved FISTA algorithm, which is significantly improved compared with the traditional method. The phase error is controlled within an acceptable range, while the effective signal energy is retained and the noise energy is greatly attenuated. (3) The present application provides efficient and high-precision reliable technical support for passive geophysical inversion, resource exploration and other fields of seismic background noise imaging technology. BRIEF DESCRIPTION OF DRAWINGS

[0016] In order to more clearly illustrate the technical solutions of the present application, the following briefly introduces the drawings needed in the embodiments. Obviously, the drawings described in the following are only some embodiments of the present application, and other drawings can also be obtained by those skilled in the art without creative labor.

[0017] Figure 1 The method flowchart of the embodiment of the present application is shown in the figure. Figure 2 The reconstructed noise cross-correlation signal waveform, front and back energy, error distribution comparison chart of the embodiment of the present application is shown in the figure, wherein, Figure 2 (a) in the figure represents the input original noise cross-correlation function schematic diagram, Figure 2 (b) in the figure represents the signal area schematic diagram, Figure 2 (c) in the figure represents the FK domain energy difference diagram, Figure 2 (d) in the figure represents the reconstructed noise cross-correlation function schematic diagram, Figure 2 (e) in the figure represents the energy difference diagram before and after reconstruction, Figure 2 (f) in the figure represents the error distribution diagram. Figure 3 FK domain spectrum of the noise cross-correlation function before and after reconstruction, the component reserved after reconstruction of the FK domain, and the filter mask diagram of the embodiment of the application, wherein, Figure 3 (a) in (a) represents the original noise cross-correlation function FK spectrum, Figure 3 (b) in (b) represents the FK spectrum after CS reconstruction, Figure 3 (c) in (c) represents the FK denoised spectrum, Figure 3 (d) in (d) represents the FK filter mask; Figure 4 The power spectrum density comparison and frequency suppression effect detection schematic diagram of the waveform at the station interval of 10.02 km, 20.54 km and 31.05 km of the original noise cross-correlation function, the compressed sensing reconstructed noise cross-correlation function and the FK filtered noise cross-correlation function; wherein, Figure 4 (a) in (a) represents the power spectrum density comparison schematic diagram of the waveform at the station interval of 10.02 km of the three noise cross-correlation functions, Figure 4 (b) in (b) represents the frequency suppression effect detection schematic diagram of the waveform at the station interval of 10.02 km of the three noise cross-correlation functions, Figure 4 (c) in (c) represents the power spectrum density comparison schematic diagram of the waveform at the station interval of 20.54 km of the three noise cross-correlation functions, Figure 4 (d) in (d) represents the frequency suppression effect detection schematic diagram of the waveform at the station interval of 20.54 km of the three noise cross-correlation functions, Figure 4 (e) in (e) represents the power spectrum density comparison schematic diagram of the waveform at the station interval of 31.05 km of the three noise cross-correlation functions, Figure 4 (f) in (f) represents the frequency suppression effect detection schematic diagram of the waveform at the station interval of 31.05 km of the three noise cross-correlation functions; Figure 5 The waveform amplitude comparison and detail window amplification schematic diagram of the original noise cross-correlation function, the compressed sensing reconstructed noise cross-correlation function and the FK filtered noise cross-correlation function at the station interval of 2.5 km, 4.01 km and 12.02 km of the embodiment of the application; wherein, Figure 5 (a) in (a) represents the waveform amplitude comparison schematic diagram of the three noise cross-correlation functions at the station interval of 2.5 km, Figure 5 (b) in (b) represents the detail window amplification schematic diagram of the three noise cross-correlation functions at the station interval of 2.5 km, Figure 5 (c) in (c) represents the waveform amplitude comparison schematic diagram of the three noise cross-correlation functions at the station interval of 4.01 km, Figure 5 (d) in (d) represents the detail window amplification schematic diagram of the three noise cross-correlation functions at the station interval of 4.01 km, Figure 5Fig. 2 (e) is a schematic diagram showing the waveform amplitude comparison of the three noise cross-correlation functions at a station interval of 12.02 km, Figure 5 Fig. 2 (f) is a schematic diagram showing the detail window amplification of the three noise cross-correlation functions at a station interval of 12.02 km; Figure 6 Fig. 2 is a schematic diagram showing the waveform amplitude comparison and power spectral density comparison of the three noise cross-correlation functions at station intervals of 2.5 km, 4.01 km and 12.02 km, wherein the three noise cross-correlation functions are the original noise cross-correlation function, the compressed sensing reconstructed noise cross-correlation function and the FK filtered noise cross-correlation function of the embodiment of the present application; Figure 6 Fig. 2 (a) is a schematic diagram showing the waveform amplitude comparison of the three noise cross-correlation functions at a station interval of 2.5 km, Figure 6 Fig. 2 (b) is a schematic diagram showing the power spectral density comparison of the three noise cross-correlation functions at a station interval of 2.5 km, Figure 6 Fig. 2 (c) is a schematic diagram showing the waveform amplitude comparison of the three noise cross-correlation functions at a station interval of 4.01 km, Figure 6 Fig. 2 (d) is a schematic diagram showing the power spectral density comparison of the three noise cross-correlation functions at a station interval of 4.01 km, Figure 6 Fig. 2 (e) is a schematic diagram showing the waveform amplitude comparison of the three noise cross-correlation functions at a station interval of 12.02 km, Figure 6 Fig. 2 (f) is a schematic diagram showing the power spectral density comparison of the three noise cross-correlation functions at a station interval of 12.02 km. DETAILED DESCRIPTION

[0018] The technical solutions in the embodiments of the present application will be apparently and completely described below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by a person of ordinary skill in the art without creative labor fall within the protection scope of the present application.

[0019] In order to make the above-mentioned objects, features and advantages of the present application more apparent and easy to understand, the present application will be further described in detail below with reference to the drawings and specific embodiments.

[0020] Embodiment One In the present embodiment, as shown in Fig. 1, a noise cross-correlation function reconstruction and purification method of seismic background noise comprises the following steps: Figure 1 Collecting original noise signals, using the sparse characteristics of the noise cross-correlation function of seismic background noise in the effective frequency band, combining the compressed sensing theory to non-uniformly downsample the original noise signals to obtain undersampling data.

[0021] ​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.

[0022] 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.

[0023] 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.

[0024] 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.

[0025] 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.

[0026] Specifically, methods for obtaining the reconstructed noise cross-correlation function include: calculating 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 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: 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 non-differentiable 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: in, x k The algorithm is represented by the first... k The sparse solution vector obtained after step iterations ensures the sparsity of the solution through a soft thresholding operation (near-end mapping); denotes a soft thresholding operator that shrinks each element of the input vector towards zero; z k denotes the k auxiliary variable (or extrapolation point) in the x k together form the iteration sequence; α denotes the step size, usually taken as and can be optimized by line search; t k denotes the sequence of time parameters that controls the Nesterov momentum acceleration in the FISTA algorithm, usually initialized as t 0=1; for high SNR, the initial value is usually set as ; the convergence threshold is set as The solution is obtained by the alternating iteration of the Nesterov accelerated gradient descent and the soft thresholding shrinkage. After 200 iterations, the final output sparse solution x200 is obtained. When the iteration reaches 200, the solution quality is usually improved effectively by increasing the iteration number: , at this time, the iteration number is usually increased to overcome the local minimum value. The 200 iterations can stably converge to the Gaussian white noise. The pulse noise needs 300+ iterations. The convergence rate of the FISTA makes the 200 iterations equivalent to 800-1000 iterations of the ISTA (Iterative Shrinkage-Thresholding Algorithm). Finally, the reconstructed noise cross-correlation function is obtained.

[0027] The reconstructed noise cross-correlation function is subjected to the frequency-wavenumber transformation, and the three-component separation is realized by using the FK filter mask to obtain the separated FK domain effective signal.

[0028] The separated FK domain effective signal includes: an effective surface wave signal, high-frequency noise, and a 0-time in-phase time delay body wave.

[0029] In this embodiment, the time-domain noise cross-correlation function is converted to the FK domain. The adaptive filter mask is constructed by using the distribution difference of the surface wave (low frequency and low wavenumber), high-frequency noise (wideband and high wavenumber), and the body wave (0-time in-phase axis, which is represented as a specific linear track) in the frequency-wavenumber space. The adaptive filter mask retains the effective surface wave through the low wavenumber passband, suppresses the high-frequency noise through the high wavenumber bandstop, and designs a notch filter for the wavenumber characteristics of the body wave interference. Finally, the separated time-domain signal is obtained through the inverse transformation, the signal-to-noise ratio is improved, and the body wave interference is eliminated, thereby providing a pure data basis for subsequent wave velocity measurement.

[0030] The separated FK domain effective signal is subjected to the inverse transformation, and the purified noise cross-correlation function is finally output.

[0031] In the embodiment, the surface wave dominant signal (retaining low frequency and low wave number components) processed by the FK filtering mask is subjected to inverse frequency-wave number transformation, and is converted from the wave number-frequency joint domain back to the time-space domain; in this process, special attention is paid to the accurate reconstruction of phase information, and the windowed inverse Fourier transform or the fast FK inversion algorithm is used to ensure the waveform integrity. The purified noise cross-correlation function finally output has three characteristics: (1) the effective surface wave dispersion characteristics are completely retained, (2) the high-frequency noise components are significantly attenuated, and (3) the 0-time phase axis body wave interference is basically eliminated, thereby providing high-quality input data meeting the requirements of the Green's function theory for subsequent crustal structure inversion.

[0032] Embodiment two In the embodiment, a noise cross-correlation function reconstruction and purification system of seismic background noise comprises a sampling module, a reconstruction module, a signal separation module and an inverse transformation module.

[0033] The sampling module is used for collecting original noise signals, and the non-uniform down-sampling of the original noise signals is performed by using the sparse characteristics of the noise cross-correlation function of seismic background noise in the effective frequency band and combining the compressed sensing theory, so as to obtain under-sampling data.

[0034] The data of the noise cross-correlation function comprises the longitude and latitude of a seismic source station, the longitude and latitude of a receiving station and the inter-station distance between the two stations.

[0035] The method for performing non-uniform down-sampling comprises: adopting Poisson disc random sampling or adaptive sampling based on signal energy distribution, and simultaneously introducing the compressed sensing theory to perform non-uniform down-sampling by using the low coherence of an observation matrix and a sparse basis, so as to obtain under-sampling data.

[0036] The reconstruction module is used for reconstructing the full-band noise cross-correlation function from the under-sampling data with high precision by using a fast iterative shrinkage threshold algorithm, so as to obtain the reconstructed noise cross-correlation function.

[0037] The method for obtaining the reconstructed noise cross-correlation function comprises: calculating the under-sampling data: wherein, y represents the under-sampling data, x represents a time sequence of the original noise cross-correlation function, represents an observation matrix, n represents noise; based on the under-sampling data, the FISTA iterative optimization is calculated: wherein, F represents a target function to be optimized, x represents a regularization parameter, λ represents a regularization parameter, Represents the sparse transformation basis; Nesterov acceleration is performed during the computation of FISTA iterative optimization: in, x k The algorithm is represented by the first... k The sparse solution vector obtained after one iteration. 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 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.

[0038] 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.

[0039] 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.

[0040] 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.

[0041] Example 3 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: 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.

[0042] 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.

[0043] The device is used for noise suppression and effective wave enhancement in seismic signal processing. The seismic data denoising eliminates linear noise (such as surface wave, multiple wave, etc.) and random noise in the seismic record through FK domain filtering, and retains the effective reflection signal. The non-uniform data standardization includes processing non-uniformly distributed data during field acquisition, and realizing FK analysis with uniform spatial interval through resampling. The velocity selective filtering is a velocity band-pass filter designed according to the difference in seismic wave propagation velocity.

[0044] The non-uniform sampling device reduces the amount of data or focuses on key signal characteristics through probabilistic sampling. Data compression reduces the amount of data by proportional down-sampling, which is suitable for optimizing seismic data storage or transmission. Energy-oriented processing prioritizes high-energy segments and suppresses low-energy noise. The algorithm comparison benchmark provides uniform sampling as a reference baseline for energy-aware sampling.

[0045] The compressed sensing reconstruction framework of the fast iterative shrinkage threshold algorithm is used for sparse recovery of one-dimensional signals and is extended to processing of two-dimensional seismic trace sets. The complete signal is recovered from non-uniformly sampled data, the data acquisition cost is reduced, and the sparsity assumption of discrete cosine transform (DCT) and the fast convergence characteristics of FISTA are used to realize efficient reconstruction.

[0046] The compressed sensing reconstruction and FK domain denoising device is used for recovering high-quality signals and suppressing noise from sparsely sampled seismic data. The complete signal is recovered from non-uniformly sampled seismic data by compressed sensing technology, noise is identified and filtered in the FK domain, and effective signals are retained. It is a one-stop solution integrating data loading, sampling, reconstruction, denoising, visualization and quality evaluation.

[0047] The seismic data visualization toolset provides full-dimensional visualization support for the seismic data processing flow. Through multi-modal chart comparison, the differences in time-space-frequency characteristics of the original data, compressed sampling reconstruction results and denoising results are compared to assist in algorithm performance evaluation and parameter optimization.

[0048] Embodiment four In this embodiment, a noise cross-correlation function reconstruction and purification method based on compressed sensing fast iterative shrinkage threshold algorithm and frequency-wave number frequency division combination includes: Step S1: Obtain the noise cross-correlation function for preprocessing to obtain the noise cross-correlation function data containing the surface wave signal; Step S2: Integrate the dispersed SAC files into a time-distance panel in a unified format, facilitate subsequent FK domain analysis or other seismic signal processing, automatically calculate the epicentral distance and sort by distance, establish a spatially ordered data structure, uniformly truncate all waveform data to the same length, ensure the dimensional matching of subsequent processing, intercept a consistent time window, and apply a cosine tapered window to suppress edge effects; Step S3: Square dimension reinforcement of high-energy area difference, energy distribution converted to probability density, ensure the rationality of irregular sampling, reduce data volume by proportional down-sampling, suitable for optimization of seismic data storage or transmission, preferentially retain high-energy segments, suppress low-energy noise, provide uniform sampling as a reference baseline for energy-aware sampling. Based on the probability distribution of waveform energy value (the higher the energy, the greater the probability of being sampled), completely random uniform sampling.

[0049] Step S4: FISTA-DCT reconstruction, where the objective function is: Where, is the signal to be reconstructed, is the observed signal containing missing values (unsampled points are set to 0), is the diagonal mask matrix ( When i points are sampled), is the DCT transform matrix, λ is the sparse regularization parameter.

[0050] Based on the input observed signal and the mask M , initialize variables: Where, denotes the Nesterov momentum auxiliary variable, is the extrapolation point (auxiliary variable) in the k th iteration, used to accelerate convergence, is a time series parameter that controls the weight of momentum, dynamically adjusts the proportion of historical information in the update, balances convergence speed and stability, initial value , step size (Lipschitz constant) is (needs to calculate or estimate the spectral norm of the matrix); Based on the FISTA main iteration, for perform gradient calculation: Where, ; Based on Nesterov accelerated update: Where, denotes a temporary variable, denotes the kextrapolation point in the next iteration; proximal mapping: element-wise soft threshold function: ; update based on inverse transform: ; momentum update: ; if the maximum number of iterations is met or the relative error condition is met, the iteration is completed: .

[0051] Step S5: Automatically detect the sparse array and switch to the conservative filtering mode. Support the standardization processing of non-uniform spatial sampling data. Add a two-dimensional cosine cone window to reduce spectral leakage, zero to the power of 2 length, FFT2 to calculate the spectrum, center, speed sector mask multiplied by the spectrum, inverse FFT to restore the time domain signal, interpolate the uniform gather back to the original non-uniform spatial position, record the key parameters such as spatial sampling rate, Nyquist wave number, and velocity resolution; Step S6: Display the time-distance profile of the original / sampling / reconstruction / denoising data side by side, compare the frequency-wavenumber energy distribution of the original / reconstruction / denoising data, extract single-channel waveforms in proportion to the distance, compare the signal patterns in the three stages, extract all channel amplitudes at fixed time points, and analyze spatial consistency.

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

[0053] In the FK domain, the main noise part is eliminated by filtering mask ( Figure 3 ), Figure 4 The power spectral density (vertical coordinate is power spectral density) of the waveform of the original noise cross-correlation function, the compressed sensing reconstructed noise cross-correlation function, and the FK filtered noise cross-correlation function at the station spacing of 10.02 km, 20.54 km, and 31.05 km is compared, and the frequency suppression effect is detected. After reconstruction, while retaining the effective signal frequency, the power of the high-frequency part is reduced.

[0054] Figure 5The 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.

[0055] 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.

[0056] Example 5 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: Data processing flow includes input and output. Key steps include data loading: reading SAC files and sorting, generating uniform length data panels. Energy-aware sampling: compressive sampling of original signals in proportion.

[0057] Compressive sensing reconstruction: using FISTA algorithm for 200 iterations to reconstruct missing data. FK denoising: filtering noise by setting velocity bands, preserving effective signals. Visualization and reporting: generating spatial-temporal domain, FK spectrum, profile map, etc. analysis charts, saving as PDF report. Core algorithms and parameters include: compressive sensing regularization parameter: lambda_reg=0.01, number of iterations: iters=200 Reconstruction effect. FK domain denoising: velocity resolution. Wave number domain mask: generated by velocity band, covering 11.5% passband.

[0058] Denoised energy: ROI area energy decreased from original 4879.575 to 901.558. Full-process visualization: provides multi-dimensional spatial-temporal and frequency domain analysis tools.

[0059] Input parameters include: --erase_roi, specify time-distance rectangular area [t0, t1, d0, d1] (seconds, kilometers), used to suppress signals in specific areas (such as direct wave interference). --edge_t, time axis soft edge transition width (seconds), avoiding Gibbs effect of rectangular filtering. --edge_d, distance axis soft edge transition width (kilometers). --in (input_dir), input directory path (must contain SAC format seismic data). --out (output_dir), output directory path (result file storage location). --ratio, compressive sampling ratio (0-1), controls data downsampling rate. --lambda (lambda_reg), compressive sensing regularization parameter, balances data fidelity and sparsity (larger values result in smoother signals). --iters, maximum number of FISTA algorithm iterations. --bands, FK filter velocity band (km / s), format as nested list, used to preserve signals in specified velocity range. --f_notch, frequency notch threshold (Hz), removes low frequency noise (such as <0.1 Hz). --k_notch, wave number 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 reproducibility of sampling process. --profile_locs, waveform comparison profile location (relative distance ratio), used for visualization analysis.

[0060] The above described embodiments are only to illustrate the preferred modes of the present application, and are not intended to limit the scope of the present application. Any modification and improvement made by those skilled in the art to the technical solutions of the present application without departing from the design spirit of the present application shall fall within the protection scope of the present application as defined by the claims.

Claims

1. A method of noise cross-correlation function reconstruction and purification of seismic ambient noise, characterized in that, The method comprises the following steps: Collecting original noise signals, using the sparse characteristics of the noise cross-correlation function of seismic background noise in the effective frequency band, and combining the theory of compressed sensing to non-uniformly down-sample the original noise signals to obtain under-sampling data; Using a fast iterative shrinkage threshold algorithm to reconstruct the full-band noise cross-correlation function from the under-sampling data with high precision to obtain the reconstructed noise cross-correlation function; Performing frequency-wavenumber transformation on the reconstructed noise cross-correlation function, and using an FK filter mask to realize three-component separation to obtain the separated FK domain effective signal; Performing inverse transformation on the separated FK domain effective signal to finally output the purified noise cross-correlation function.

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

3. The method of claim 1, wherein, The method for performing the non-uniform down-sampling comprises: Using Poisson disc random sampling or adaptive sampling based on signal energy distribution, and simultaneously introducing the theory of compressed sensing, non-uniform down-sampling is performed through the low coherence of an observation matrix and a sparse basis to obtain under-sampling data.

4. The method of claim 1, wherein, The method for obtaining the reconstructed noise cross-correlation function comprises: Calculating the under-sampling 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 under-sampling data, FISTA iterative optimization is calculated: wherein F x represents an objective function to be optimized, During the calculation of the FISTA iterative optimization, Nesterov acceleration is performed: represents a regularization parameter, represents a sparse transform basis;​ Through the alternative iteration of Nesterov accelerated gradient descent and soft threshold shrinkage, the reconstructed noise cross-correlation function is obtained. wherein, x k denotes the algorithm, k denotes the sparse solution vector obtained after the iteration step, denotes the soft thresholding operator, z k denotes the iteration, k denotes the auxiliary variable in the iteration, α denotes the step size, t k denotes the sequence of time parameters controlling the Nesterov momentum acceleration in the FISTA algorithm; The separated FK domain effective signal comprises an effective surface wave signal, high-frequency noise, and a 0-time in-phase time delay body wave.

5. The method of claim 1, wherein, It comprises:

6. A system for reconstruction and purification of noise cross-correlation functions of seismic ambient noise, said system applying the method according to any one of claims 1 to 5, characterized in that, a sampling module, a reconstruction module, a signal separation module, and an inverse transformation module; The sampling module is used to collect original noise signals, use the sparse characteristics of the noise cross-correlation function of seismic background noise in the effective frequency band, and combine the theory of compressed sensing to non-uniformly down-sample the original noise signals to obtain under-sampling data; The reconstruction module is used to use a fast iterative shrinkage threshold algorithm to reconstruct the full-band noise cross-correlation function from the under-sampling data with high precision to 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 realize three-component separation to obtain the separated FK domain effective signal; The inverse transformation module is used to perform inverse transformation on the separated FK domain effective signal to finally output the purified noise cross-correlation function. The data of the noise cross-correlation function comprises the longitude and latitude of a seismic source station, the longitude and latitude of a receiving station, and the inter-station distance between the two stations.

7. The system for reconstruction and purification of noise cross-correlation function of seismic ambient noise according to claim 6, wherein, The method for performing the non-uniform down-sampling comprises:

8. The system for reconstruction and purification of noise cross-correlation function of seismic ambient noise according to claim 6, wherein, Using Poisson disc random sampling or adaptive sampling based on signal energy distribution, and simultaneously introducing the theory of compressed sensing, non-uniform down-sampling is performed through the low coherence of an observation matrix and a sparse basis to obtain under-sampling data. The method for obtaining the reconstructed noise cross-correlation function comprises:

9. The system for reconstruction and purification of noise cross-correlation function of seismic ambient noise according to claim 6, characterized in that, Calculating the under-sampling data: Based on the under-sampling data, FISTA iterative optimization is calculated: 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; During the calculation of the FISTA iterative optimization, Nesterov acceleration is performed: wherein F x represents the objective function to be optimized, ​ represents a regularization parameter, represents a sparse transform basis;​ ​ wherein, x k denotes the algorithm, k denotes the sparse solution vector obtained after the iteration step, denotes the soft thresholding operator, z k denotes the auxiliary variable in the iteration, k denotes the step size, α denotes the step size, t k denotes the sequence of time parameters controlling the Nesterov momentum acceleration in the FISTA algorithm; The reconstructed noise cross-correlation function is solved by Nesterov accelerated gradient descent and soft threshold shrinkage iterative alternation.

10. The system for reconstruction and purification of noise cross-correlation function of seismic ambient noise according to claim 6, wherein, The separated FK domain effective signal comprises an effective surface wave signal, high frequency noise and a 0-time in-phase time delay body wave.

Images