Seismic signal multi-scale and OVT denoising processing system based on time-frequency analysis

Abstract

The invention relates to the technical field of seismic signal exploration, in particular to a seismic signal multi-scale and OVT denoising processing system based on time-frequency analysis, which comprises a data acquisition module used for acquiring a seismic signal SEQ; the system is further used for collecting tide data, ocean current data, coordinates (xHF, yHF and zHF) of seismic signal collection positions, ocean depth DOcean, collection frequency fEQ, water temperature WT, salinity SAL and water pressure WP. According to the invention, the tide and ocean current influence model is established, meteorological data and hydrological data are combined, the tide water change and the ocean flow velocity fluctuation are effectively tracked, and the environmental change in the seismic signal propagation process can be timely responded; through setting a machine learning model, parameters of seismic signal analysis can be dynamically adjusted (that is, influence of tides and ocean currents is compensated), so that influence factors of seismic signals in a propagation process are considered, accuracy and reliability of data processing are improved, and accurate extraction of signal features is ensured.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic signal surveying, and in particular to a seismic signal multi-scale and OVT denoising processing system based on time-frequency analysis. Background Art

[0002] Seismic signal processing is an important direction in seismological research, and its core task is to extract and analyze seismic signals to obtain effective information. The analysis and processing of seismic signals are extremely important in earthquake monitoring, early warning, and subsequent disaster assessment. Traditional seismic signal processing methods are mostly focused on time domain analysis. Although they can effectively extract signal features under certain circumstances, the effect of time domain analysis is often limited when faced with complex noise interference and non-stationary signals.

[0003] After searching, Chinese patent number CN201910950889.2 discloses a pre-stack seismic data processing method in the OVT domain, which includes: sorting the original seismic data according to different azimuths and offset distances to obtain OVT domain gathers; decomposing the data into instantaneous amplitude and instantaneous phase in the OVT domain; selecting and processing the central surface element, and superimposing the gathers within a certain spatial range of the central surface element based on the superposition principle; repeating the above processing steps for all central surface elements to obtain the results of all processed central surface elements, and reconstructing the output phase data and amplitude data.

[0004] In the above scheme, although the original data are sorted into OVT gathers with different offset distances and different azimuth angles, the correlation between data is enhanced because the OVT gathers have the same offset distance and azimuth angle information, the data are decomposed into instantaneous amplitude and instantaneous phase in the OVT domain, and weak signal enhancement is achieved in the complex domain. The data after weak signal amplitude enhancement and phase enhancement are reconstructed, and the effective weak reflection signal masked by noise is enhanced on the basis of relative data fidelity, thereby improving the signal-to-noise ratio of pre-stack data as a whole, so that the effects of pre-stack gathers, stacking and imaging profiles are significantly improved.

[0005] However, in actual applications, conventional seismic signal processing technology often faces serious challenges in complex geological environments, such as areas with faults, folds or salt domes. The complex geological structures in these areas lead to the diversity of seismic wave propagation characteristics, which in turn increases the interference between signals and noise. In addition, in the marine environment, changes in tides and ocean currents will also affect the propagation speed and propagation path of seismic waves. This dynamic change causes the time-frequency characteristics of the signal to change continuously, increasing the difficulty of processing.

[0006] Therefore, a seismic signal multi-scale and OVT denoising processing system based on time-frequency analysis is proposed to solve the above problems. Summary of the invention

[0007] Technical issues solved

[0008] In view of the above-mentioned shortcomings of the prior art, the present invention provides a multi-scale and OVT denoising processing system for seismic signals based on time-frequency analysis, which can effectively solve the problem of considering the influence of geological environment characteristics when processing seismic data based on the OVT domain in the prior art, especially the interference problem of tidal and current changes in the marine environment.

[0009] Technical Solution

[0010] To achieve the above objectives, the present invention is implemented through the following technical solutions:

[0011] The present invention provides a multi-scale and OVT denoising processing system for seismic signals based on time-frequency analysis, comprising a data acquisition module for acquiring seismic signals S EQ ; It is also used to collect tidal data and ocean current data, as well as the coordinates of the location where the seismic signal is collected (x HF ,y HF ,z HF ), Ocean Depth D O cean, acquisition frequency f EQ , water temperature WT, salinity SAL and water pressure WP;

[0012] The data extraction module is used to perform wavelet transform on the seismic signal to obtain the frequency component, extract the high-frequency component in the frequency component for processing, construct the processed high-frequency component and low-frequency component into an effective signal, and calculate the signal-to-noise ratio (SNR) of the effective signal;

[0013] The data separation module is used to separate the independent components IC from the effective signal and merge all the independent components IC into the reconstructed signal S re , and calculate the reconstruction error E re ;

[0014] The data processing module builds a tide and current impact model based on the reconstructed signal to obtain a correction signal;

[0015] Processing enhancement module, based on the coordinates of the seismic signal acquisition location (x HF ,y HF ,z HF ), Ocean Depth D Ocean , acquisition frequency f EQ , water temperature WT, salinity SAL and water pressure WP are used to update the tide and current impact model and obtain the enhanced signal S Reinf (T(D Ocean ,f EQ )).

[0016] Furthermore, for the seismic signal S EQThe way to perform wavelet transform is:

[0017] According to the preset decomposition scale number n, the seismic signal is decomposed into n frequency components F(S EQ ), the decomposition formula is: In the formula, S j,k is the wavelet coefficient, which represents the local features at a given scale j and position k; S EQ (t) is the seismic signal at time t; is the wavelet function; a is the scale parameter; k is the displacement parameter; t is the time;

[0018] Define the frequency component F(S EQ ) is calculated as: Where J represents the wavelet transform of the seismic signal S EQ The maximum scale for decomposition; K represents the wavelet coefficient S at each scale j j,k The number of locations; F(j) represents the frequency value corresponding to a specific scale j, Where C is a constant, a j is the scale factor associated with scale j.

[0019] Furthermore, the high frequency components are processed in the following manner:

[0020] Extract the frequency component F(S EQ ) high-frequency components, the extraction formula is:

[0021] In the formula, H(F(S EQ )) represents the transfer function of the high-pass filter; F C Preset frequency value for high pass filter;.

[0022] Calculate the mean and standard deviation of the high-frequency components using the following formula:

[0023]

[0024] Where μ is the mean value of the high-frequency component; N is the total number of signal values ​​in the high-frequency component; is the signal value of the ith high-frequency component; σ is the standard deviation of the high-frequency component;

[0025] Set the dynamic artifact threshold TH ART =μ+κ·σ; where κ is a constant coefficient;

[0026] The mean and standard deviation of the high-frequency components are updated based on the sliding window. The updating method is:

[0027] Set the size of the sliding window and the initial mean μ init , standard deviation σ initand the corresponding dynamic artifact threshold TH ARTinit Slide the sliding window forward by 1 sample, remove the earliest sample in the sliding window, add the next sample, calculate the new mean and standard deviation, and update the dynamic artifact threshold to get the mean μ of the next sliding window init+1 , standard deviation σ init+1 and the dynamic artifact threshold TH ARTinit+1 ; Repeat the above process until all samples are processed; finally get the processed high frequency component COMP h ' igh .

[0028] Furthermore, the calculation formula of the signal-to-noise ratio SNR of the effective signal is:

[0029] Where P signal is the signal power of the effective signal, and the calculation formula is: Where Q is the number of samples of the valid signal, q is the serial index of the valid signal; P noise is the noise power, and the calculation formula is: Where P is the number of noise samples, Noise p is the remaining noise sample, obtained based on the above waveform comparison and analysis, and p is the serial index of the noise.

[0030] Furthermore, the method of separating out the independent components includes:

[0031] The signal-to-noise ratio SNR is greater than the signal-to-noise ratio threshold TH SNR The effective signal S eff Extract, center and whiten it to get the effective signal after processing After processing, the effective signal Separate and extract the independent component IC by:

[0032] Set the initial weight matrix W to convert the processed valid signal Mapped to the matrix of the independent component space; Update the initial weight matrix W until it converges; the update formula is: Where W new is the new weight matrix obtained after the current iteration, W old is the old weight matrix obtained in the previous iteration; g(...) is the nonlinear activation function; g′(...) is the derivative of the nonlinear activation function; ψ{...} is the expectation operator;

[0033] The formula for extracting independent components is: Where W final is the final weight matrix after iterative update.

[0034] Furthermore, the calculation of the reconstruction error E re The formula is:

[0035] Where M is the effective signal after processing The number of samples.

[0036] Furthermore, the method of constructing the tide and current impact model is:

[0037] Based on tidal data S Tide and ocean current data S Ocea Calculate the change in signal propagation speed ΔV = a·S Tide +b·S Ocean +c; where ΔV represents the change in the propagation velocity of seismic signals under the influence of tides and ocean currents; a represents the influence of tidal changes on the propagation velocity of seismic signals, b represents the influence of ocean currents on the propagation velocity of seismic signals, and c represents the basic propagation velocity;

[0038] Calculate the adjusted propagation speed V based on the signal propagation speed change ΔV adj =ΔV+V init , where V init is the initial propagation speed;

[0039] Based on the adjusted propagation speed V adj For the reconstructed signal S re The wave equation is updated to obtain Where P EQ represents the pressure fluctuation of the reconstructed signal; is the Laplace operator; the pressure fluctuation P of the reconstructed signal EQ Perform Hilbert transform to obtain the complex signal of the reconstructed signal Where P EQ (t) is the pressure fluctuation at time point t; P EQ The Hilbert transform result of ; χ is an imaginary unit; the instantaneous amplitude A(t) and instantaneous phase φ(t) are extracted from the complex signal X(t) of the reconstructed signal, and the extraction formula is: A(t) = |X(t), φ(t) = arg(X(t));

[0040] The correction signal is calculated using the instantaneous amplitude A(t) and instantaneous phase φ(t). The calculation formula is: S Corr (t) = A(t)·e jφ(t) ; In the formula, e jφ(t) is the complex exponential form; S Corr (t) is the correction signal.

[0041] Furthermore, the training method of the tide and current impact model is as follows:

[0042] Define a recursive neural network model as the basic structure of the tide and current impact model, the basic structure includes an input layer, a hidden layer, and an output layer;

[0043] Introduce the parameter vector time step d, and transform the vector time step d and tidal data S Tide and ocean current data S Ocean Combined, we get the fused feature vector X d ={(S Tide ,S Ocean ) d-1 ,(S Tide ,S Ocean ) d-2 ,...,(S Tide ,S Ocean ) d-s}, where s is the number of vector time steps d;

[0044] Collect historical tide data S within a fixed period of time in the past Tide and historical ocean current data S Ocean , and correspondingly construct the historical fusion feature vector Fusion of history into feature vector Expand by vector time steps to form sequence data as training input for the tide and current impact model;

[0045] Initialize the network parameters of the tide and current impact model; the network parameters include the weight matrix W from the input layer to the hidden layer x , the recurrent weight matrix W of the hidden layer h , the bias vector b of the hidden layer h , the weight matrix W from the hidden layer to the output layer o and the bias vector b of the output layer o ;

[0046] Define the loss function for the model affected by tides and currents

[0047] In the formula, S is the length of the sequence data; Y true,d is the true value of the dth vector time step; Y pred,d is the predicted value of the dth vector time step;

[0048] The sequence data is passed to the input layer, and then passes through the hidden layer and the output layer in turn, the predicted value is output, the corresponding loss function value is calculated, and the error gradient is back-propagated from the output layer to the input layer; according to the error gradient, the network parameters are updated using the optimization algorithm; the sequence data is repeatedly trained until the tide and current impact model converges or reaches the preset number of iterations, which means that the tide and current impact model training is completed.

[0049] Furthermore, the input layer is used to receive the fused feature vector X d As input; the hidden layer consists of several state vectors, and the calculation formula for each state vector is: d =f(W h ·h d-1 +W x ·X d +b h );where h d is the hidden state of the current vector at time step d, and f(...) is the activation function;

[0050] The output layer calculates the output results of the tide and current impact model based on the output of the hidden layer. The calculation formula is: S Corr (t) = η(W O ·h d +b o ), where η(...) is the activation function of the output layer.

[0051] Furthermore, the method for updating the tide and current impact model is as follows:

[0052] Calculate the ocean depth D at the current acquisition location Ocean The sound wave propagation speed V(D Ocean ,f EQ )=V0+ΔV(D Ocean );where V0 is the speed of sound at unit speed; ΔV(D Ocean ) is the ocean depth D at the current acquisition location Ocean The change in the speed of sound waves under water temperature, salinity and pressure is calculated as: ΔV(D Ocean )=μ1·WT+μ2·SAL+μ3·WP; where μ1 is the influence coefficient of water temperature on the sound wave velocity, μ2 is the influence coefficient of salinity on the sound wave velocity, and μ3 is the influence coefficient of water pressure on the sound wave velocity;

[0053] Calculate signal propagation time Where, the signal propagation time T Prop (D Ocean ,f EQ ) represents the signal propagation time from the source to the receiver; Dist represents the distance from the source to the receiver, Where (x RX ,y RX ,z RX ) are the coordinates of the receiver;

[0054] The change in the speed of the sound wave ΔV(D Ocean ) and signal propagation time T Prop (D Ocean ,f EQ) is returned to the tide and current influence model to correct the signal S Corr (t) is updated to obtain the enhanced signal S Reinf (T(D Ocean ,f EQ )).

[0055] Beneficial Effects

[0056] Compared with the prior art, the technical solution provided by the present invention has the following beneficial effects:

[0057] The present invention transforms the seismic signal S EQ After decomposing into signal components of multiple scales (high-frequency components and low-frequency components), the high-frequency part (such as instantaneous changes in seismic signals) and the low-frequency part (such as the overall trend of events) can be processed in a targeted manner, making the subsequent analysis of seismic signals more accurate; and by extracting high-frequency components, the background noise in the seismic signal can be effectively removed, ensuring that subsequent analysis focuses only on key information, thereby improving the quality of the data;

[0058] This scheme receives the effective signal S eff Separation can effectively remove the interference caused by multipath propagation, thereby improving the effective signal S eff The clarity of the signal can be effectively identified and extracted. re , reducing the interference and improving the clarity and interpretability of the signal; this means that in subsequent analysis, the reconstructed signal S re The main features of the seismic signal are more obvious, which helps to accurately identify the interference of potential geological structures on seismic signals;

[0059] This solution establishes a tidal and current impact model, combines meteorological data with hydrological data, effectively tracks tidal changes and ocean velocity fluctuations, and can respond promptly to environmental changes during the propagation of seismic signals. Through the setting of machine learning models, the parameters of seismic signal analysis can be dynamically adjusted (i.e., by compensating for the effects of tides and currents), so that the influencing factors of seismic signals in the propagation process can be taken into account, thereby improving the accuracy and reliability of data processing and ensuring the accurate extraction of signal features. BRIEF DESCRIPTION OF THE DRAWINGS

[0060] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the prior art descriptions are briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention, and for ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0061] Figure 1Schematic diagram of the system structure in an embodiment of the present invention. DETAILED DESCRIPTION

[0062] In order to make the purpose, technical solution and advantages of the embodiments of the present invention clearer, the technical solution in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0063] The present invention will be further described below in conjunction with the embodiments.

[0064] Embodiment 1:

[0065] See attached Figure 1 This case proposes a multi-scale and OVT denoising processing system for seismic signals based on time-frequency analysis, including a data acquisition module, a data extraction module, a data separation module, a data processing module and a data enhancement module; among which:

[0066] The data acquisition module acquires seismic signals based on seismic sensors arranged in the marine exploration area. Its sampling frequency is usually set to 1000Hz to 400Hz to ensure that sufficient seismic signals are captured. It is also used to collect tidal data and ocean current data; as well as the coordinates of the seismic signal acquisition location, ocean depth, acquisition frequency, water temperature, salinity and water pressure;

[0067] The data extraction module is used to perform wavelet transform on the seismic signal, decompose the seismic signal into n frequency components, extract the high-frequency components from the frequency components, detrend and remove artifacts from the high-frequency components, remove the DC components (i.e., the average value of the seismic signal), and identify and remove high-frequency components that exceed the normal range; construct the processed high-frequency components and low-frequency components into effective signals, and calculate their signal-to-noise ratio to verify the quality of the effective signals;

[0068] The data separation module is used to decompose the effective signal to obtain the processed effective signal, then separate different independent components from the processed effective signal, and merge all independent components into a reconstructed signal to reduce the interference caused by multipath propagation. Finally, the reconstruction error E is calculated. re Evaluate the reconstructed signal S re quality;

[0069] The data processing module builds a tide and current impact model based on the reconstructed signal to obtain a correction signal;

[0070] The processing enhancement module updates the tide and current impact model based on the coordinates of the seismic signal acquisition location, ocean depth, acquisition frequency, water temperature, salinity and water pressure to obtain the enhanced signal S Reinf (T(D Ocean ,f EQ )).

[0071] Specifically, in a complex marine geological environment (such as a geological environment with faults or salt domes), seismic signals may be affected by strong background noise, resulting in the overlap of background noise and effective signals during detrending and artifact removal, making it difficult to accurately identify and remove DC components during detrending; artifacts may appear as instantaneous high-amplitude changes, and the set threshold may not be able to effectively distinguish between real signals and artifacts; in this scheme, the seismic signal is preprocessed as follows:

[0072] Based on wavelet transform, the seismic signal S EQ Decompose the seismic signal into n frequency components F(S EQ ), the decomposition formula is:

[0073] In the formula, S j,k is the wavelet coefficient, which represents the local characteristics at a given scale j and position k and can be used to analyze the time and frequency characteristics of seismic signals; S EQ (t) is the seismic signal at time t; is a wavelet function, which determines the characteristics of the transformation and the way to process seismic signals. Daubechies wavelet or Morlet wavelet can be used; a is a scale parameter, which is used to control the width and expansion of the wavelet; k is a displacement parameter, which represents the wavelet function Position in the time domain; t is the time instant; is the normalization factor used to normalize the wavelet coefficient S j,k , to ensure good comparability between coefficients at different scales;

[0074] Define the frequency component F(S EQ ) is the representation of the seismic signal in the frequency domain, and its calculation formula is:

[0075] Where J represents the wavelet transform of the seismic signal S EQ The maximum scale for decomposition determines how many scales of frequency components can be extracted, affecting the resolution between high and low frequencies; K represents the wavelet coefficient S at each scale j. j,k The number of locations reflects the time localization ability of the signal; F(j) represents the frequency value corresponding to a specific scale j, Where C is a constant, a jis the scale factor associated with scale j;

[0076] Extract the frequency component F(S EQ ) and retain the low-frequency components, which helps to separate the effective information and noise in the seismic signal. The extraction formula is: In the formula, H(F(S EQ )) represents the transfer function of the high-pass filter, which is used to determine how the seismic signal with high-frequency components passes through the filter; F C The frequency value preset for the high-pass filter is used to distinguish the high-frequency component COMP high and low frequency components COMP low ;

[0077] In summary, the seismic signal S is transformed into EQ After decomposing into signal components of multiple scales (high-frequency components and low-frequency components), the high-frequency part (such as instantaneous changes in seismic signals) and the low-frequency part (such as the overall trend of events) can be processed in a targeted manner, making the subsequent analysis of seismic signals more accurate; and by extracting high-frequency components, the background noise in the seismic signal can be effectively removed, ensuring that subsequent analysis only focuses on key information, thereby improving the quality of the data.

[0078] Set the dynamic artifact threshold to remove noise and interference in the high-frequency component and obtain the processed high-frequency component COMP h ' igh The setting method is;

[0079] Calculate the mean and standard deviation of the high-frequency components using the following formula:

[0080]

[0081] Where μ is the mean value of the high-frequency component; N is the total number of signal values ​​in the high-frequency component; is the signal value of the ith high-frequency component; σ is the standard deviation of the high-frequency component;

[0082] Set the dynamic artifact threshold based on the mean and standard deviation of the high-frequency components: TH ART =μ+κ·σ; where TH ART is the dynamic artifact threshold; κ is the constant coefficient;

[0083] The mean and standard deviation of the high-frequency components are updated based on the sliding window. The updating method is:

[0084] Set the size of the sliding window and the initial mean μ init , standard deviation σ init and the corresponding dynamic artifact threshold TH ARTinitSlide the sliding window forward by 1 sample, remove the earliest sample in the sliding window, add the next sample, calculate the new mean and standard deviation, and update the dynamic artifact threshold to get the mean μ of the next sliding window init+1 , standard deviation σ init+1 and the dynamic artifact threshold TH ARTinit+1 ; Repeat the above process until all samples are processed; finally get the processed high frequency component COMP h ' igh ;

[0085] The dynamic artifact threshold TH is updated in real time through a sliding window ART , can automatically adapt to the characteristic changes of high-frequency components to ensure the set dynamic artifact threshold TH ART Always reflects the temporary high-frequency component status, thereby improving the accuracy and effect of artifact removal;

[0086] The processed high frequency component COMP h ' igh and low frequency components COMP low Constructed into a valid signal S eff ; Compare the seismic signal waveform before processing with the effective signal waveform, identify the artifact areas to be removed, observe the relationship between these artifact areas and the effective signal, ensure that the removed part is noise rather than effective signal, and ensure that the effective signal is not mistakenly removed;

[0087] Calculate the signal-to-noise ratio (SNR) of the effective signal using the following formula:

[0088] Where P signal is the signal power of the effective signal, and the calculation formula is: Where Q is the number of samples of the valid signal, q is the serial index of the valid signal; P noise is the noise power, and the calculation formula is: Where P is the number of noise samples, Noise p is the remaining noise sample, obtained based on the above waveform comparison analysis, and p is the serial index of the noise;

[0089] Based on the signal-to-noise ratio SNR, the seismic signal S EQ The processing results (i.e., effective signals) are evaluated. A higher signal-to-noise ratio indicates that the intensity of the effective signal is significantly higher than the noise, ensuring the signal quality, thereby reducing the risk of misjudgment in subsequent analysis and providing feedback basis for adjusting and optimizing subsequent processing steps.

[0090] Furthermore, in this solution, the method of extracting and decomposing the effective signal and separating different independent components therefrom is as follows:

[0091] The signal-to-noise ratio SNR is greater than the preset signal-to-noise ratio threshold TH SNR The effective signal S eff Extract, center and whiten it to get the effective signal after processing The centering process refers to removing the mean of the effective signal so that the mean value of the effective signal is zero. This step helps to eliminate the DC component in the signal and facilitates subsequent analysis. The whitening process refers to making the variances of different effective signals equal and the effective signals uncorrelated.

[0092] Based on FastlCA algorithm, the effective signal after processing Separate to obtain independent components; the separation method is:

[0093] Set the initial weight matrix W to convert the processed valid signal The matrix mapped to the independent component space; the initial weight matrix W is updated until it converges (that is, the difference between the new weight matrix and the old weight matrix is ​​less than the preset threshold, thereby achieving convergence); the update formula is: Where W new is the new weight matrix obtained after the current iteration, W old is the old weight matrix obtained in the previous iteration; g(...) is a nonlinear activation function, such as a hyperbolic tangent or logistic sigmoid function; g′(...) is the derivative of the nonlinear activation function, which is used to calculate the rate of change of the activation function to the processed effective signal; ψ{...} is the expectation operator, which is used to calculate the processed effective signal of the input The average result after activation function processing; the formula for extracting independent components is: Where IC is the independent component and is the effective signal after processing. The signal matrix separated from final is the final weight matrix after iterative update, indicating the effective signal after processing To the mapping relationship of independent components; merge all independent components IC to obtain the reconstructed signal S re ; and calculate the reconstruction error E re , the calculation formula is: Where M is the effective signal after processing The number of samples; according to the reconstruction error E re The size of the reconstructed signal S is evaluated re quality.

[0094] In summary, for the received valid signal S eff Separation can effectively remove the interference caused by multipath propagation, thereby improving the effective signal S effThe clarity of the signal can be effectively identified and extracted. re , reducing the interference and improving the clarity and interpretability of the signal; this means that in subsequent analysis, the reconstructed signal S re The main features of the seismic signal are more obvious, which helps to accurately identify the interference of potential geological structures on seismic signals.

[0095] Specifically, the method of constructing the tide and current impact model is as follows:

[0096] Based on tidal data S Tide and ocean current data S Ocean Calculate the change in signal propagation speed ΔV, the calculation formula is: ΔV = a·S Tide +b·S Ocean +c; where ΔV represents the change in the propagation velocity of the seismic signal under the influence of tides and ocean currents; a, b, and c are all constants, a represents the influence of tidal changes on the propagation velocity of the seismic signal, b represents the influence of ocean currents on the propagation velocity of the seismic signal, and c represents the basic propagation velocity; under hydrological and meteorological conditions, the influence of tide level and ocean current velocity on the propagation velocity can be reasonably approximated as a linear relationship, and the above formula can provide sufficient accuracy to calculate the signal propagation velocity change ΔV without introducing unnecessary complexity;

[0097] Calculate the adjusted propagation speed V based on the signal propagation speed change ΔV adj , the calculation formula is: V adj =ΔV+V init , where V init is the initial propagation speed, which is usually obtained based on environmental data or experimental data;

[0098] Based on the adjusted propagation speed V adj For the reconstructed signal S re The wave equation is updated to obtain Where P EQ The pressure fluctuations representing the reconstructed signal reflect the changes in pressure in the medium during the propagation of seismic waves; is the Laplace operator, which represents the change of pressure fluctuation of the reconstructed signal in space; The second-order derivative of the pressure fluctuation of the reconstructed signal with respect to time t describes the rate of change of the pressure fluctuation over time;

[0099] The pressure fluctuation P of the reconstructed signal EQ Perform Hilbert transform to obtain the complex signal of the reconstructed signal Where P EQ (t) is the pressure fluctuation at time point t; P EQThe Hilbert transform result provides the phase information of the pressure fluctuation; χ is the imaginary unit;

[0100] Extract the instantaneous amplitude A(t) and instantaneous phase φ(t) from the complex signal X(t) of the reconstructed signal. The extraction formula is: A(t) = |X(t)|, φ(t) = arg(X(t)); the instantaneous amplitude A(t) reflects the strength or energy of the reconstructed signal and can reveal the attenuation and enhancement of the seismic wave; the instantaneous phase φ(t) represents the phase position of the reconstructed signal, which is affected by the adjusted propagation velocity in the wave equation and can reveal the propagation characteristics and phase relationship of the reconstructed signal.

[0101] The correction signal is calculated using the instantaneous amplitude A(t) and instantaneous phase φ(t). The calculation formula is: S Corr (t) = A(t)·e jφ(t) ; In the formula, e jφ(t) It is a complex exponential form, representing the phase information; S Corr (t) is the correction signal, which is composed of instantaneous amplitude and phase, reflecting the real fluctuation of the seismic signal in the time domain.

[0102] Define a recursive neural network model as the basic structure of the tide and current impact model, the basic structure includes an input layer, a hidden layer, and an output layer;

[0103] Introduce the parameter vector time step d, and transform the vector time step d and tidal data S Tide and ocean current data S Ocean Combined, we get the fused feature vector X d ={(S Tide ,S Ocean ) d-1 ,(S Tide ,S Ocean ) d-2 ,...,(S Tide ,S Ocean ) d-s}, where s is the number of vector time steps d;

[0104] The input layer is used to receive the fused feature vector X d As input; the hidden layer consists of several state vectors, and the calculation formula for each state vector is: d =f(W h ·h d-1 +W x ·X d +b h );where h d is the hidden state of the current vector at time step d, f(...) is the activation function, W x is the weight matrix from the input layer to the hidden layer, W his the recurrent weight matrix of the hidden layer, b h is the bias vector of the hidden layer;

[0105] The output layer is based on the output of the hidden layer, i.e., the hidden state h at the current vector time step d. d , the output of the tide and current impact model is calculated, which is the correction signal S Corr (t), the calculation formula is: S Corr (t) = η(W O ·h d +b o ), where η(...) is the activation function of the output layer, W o is the weight matrix from the hidden layer to the output layer, b o is the bias vector of the output layer;

[0106] Collect historical tide data S within a fixed period of time in the past Tide and historical ocean current data S Ocean , and the corresponding fusion feature vector is constructed, recorded as the historical fusion feature vector Fusion of history into feature vector Expand by vector time steps to form sequence data as training input for the tide and current impact model;

[0107] Initialize the network parameters of the tide and current impact model; the network parameters include the weight matrix W from the input layer to the hidden layer x , the recurrent weight matrix W of the hidden layer h , the bias vector b of the hidden layer h , the weight matrix W from the hidden layer to the output layer o and the bias vector b of the output layer o ;

[0108] Define the loss function for the model affected by tides and currents

[0109] In the formula, S is the length of the sequence data; Y true,d is the true value of the dth vector time step, that is, the historical correction signal S′ Corr (t); Y pred,d is the predicted value of the dth vector time step, that is, the predicted correction signal S Corr (t);

[0110] The sequence data is passed to the input layer, and then passes through the hidden layer and the output layer in turn to output the predicted correction signal S Corr(t), and calculate the corresponding loss function value, and back-propagate the error gradient along the output layer to the input layer; according to the error gradient, use the optimization algorithm (such as gradient descent, Adam, etc.) to update the network parameters to reduce the value of the loss function; repeat the training of the sequence data until the tide and current impact model converges (that is, the value of the loss function no longer changes) or reaches the preset number of iterations, which means that the tide and current impact model training is completed.

[0111] In summary, by establishing a tidal and ocean current impact model and combining meteorological data with hydrological data, we can effectively track tidal changes and ocean current velocity fluctuations, and make timely responses to environmental changes during the propagation of seismic signals. And through the setting of machine learning models, we can dynamically adjust the parameters of seismic signal analysis (that is, by compensating for the effects of tides and ocean currents), so that the influencing factors of seismic signals in the propagation process can be taken into account, thereby improving the accuracy and reliability of data processing and ensuring the accurate extraction of signal features.

[0112] It is worth noting that in this case, when the seismic signal is collected, the different levels of seawater depth and frequency response will affect the seismic signal analysis, resulting in uncertainty in signal processing; since the seawater depth is related to the speed of sound wave propagation, deeper or gradual water layers will affect the propagation and reception performance of the source signal at different frequencies. Therefore, by optimizing the processing capabilities of seismic signals at different seawater depths and frequency responses, the uncertainty in the signal analysis process can be reduced and the clarity and accuracy of the signal can be improved; specifically:

[0113] When the data acquisition module collects seismic signals, it also marks the coordinates of the seismic signal acquisition location (x HF ,y HF ,z HF ), Ocean Depth D O cean, acquisition frequency f EQ , water temperature WT, salinity SAL and water pressure WP; calculate the ocean depth D at the current acquisition location Ocean The sound wave propagation speed V(D Ocean ,f EQ )=V0+ΔV(D Ocean );where V0 is the speed of sound at unit speed; ΔV(D Ocean ) is the ocean depth D at the current acquisition location Ocean The change in the speed of sound waves under water temperature, salinity and pressure is calculated as: ΔV(D Ocean )=μ1·WT+μ2·SAL+μ3·WP; where μ1 is the influence coefficient of water temperature on the sound wave velocity, μ2 is the influence coefficient of salinity on the sound wave velocity, and μ3 is the influence coefficient of water pressure on the sound wave velocity. μ1, μ2 and μ3 are all obtained through experimental data;

[0114] Calculate signal propagation time Where, the signal propagation time T Prop (D Ocean ,f EQ ) represents the signal propagation time from the source to the receiver; Dist represents the distance from the source to the receiver, Where (x RX ,y RX ,z RX ) are the coordinates of the receiver;

[0115] The change in the speed of the sound wave ΔV(D Ocean ) and signal propagation time T Prop (D Ocean ,f EQ ) is returned to the tide and current influence model to correct the signal S Corr (t) is updated to obtain the enhanced signal S Reinf (T(D Ocean ,f EQ ));

[0116] By considering the influence of the coordinates of the seismic signal acquisition location, ocean depth, acquisition frequency, water temperature, salinity and water pressure on the transmission of sound waves, the processing of seismic signals can be made more accurate, ensuring the adaptability of the system to the processing of seismic signals collected in different environments, and significantly improving the accuracy and effectiveness of seismic signal analysis.

[0117] The above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements will not cause the essence of the corresponding technical solutions to deviate from the protection scope of the technical solutions of the embodiments of the present invention.

Claims

1. A multi-scale and OVT denoising processing system for seismic signals based on time-frequency analysis, characterized in that: Including data acquisition module for collecting seismic signals S EQ ; It is also used to collect tidal data and ocean current data, as well as the coordinates of the location where the seismic signal is collected (x HF ,y HF ,z HF ), Ocean Depth D O cean, acquisition frequency f EQ , water temperature WT, salinity SAL and water pressure WP; The data extraction module is used to perform wavelet transform on the seismic signal to obtain the frequency component, extract the high-frequency component in the frequency component for processing, construct the processed high-frequency component and low-frequency component into an effective signal, and calculate the signal-to-noise ratio (SNR) of the effective signal; Data separation module, used to separate independent component IC from valid signal, and merge all independent component IC into S re , and calculate the reconstruction error E re ; The data processing module builds a tide and current impact model based on the reconstructed signal to obtain a correction signal; Processing enhancement module, based on the coordinates of the seismic signal acquisition location (x HF ,y HF ,z HF ), Ocean Depth D Ocean , acquisition frequency f EQ , water temperature WT, salinity SAL and water pressure WP are used to update the tide and current impact model and obtain the enhanced signal S Reinf (T(D Ocean ,f EQ )).

2. The seismic signal multi-scale and OVT denoising processing system based on time-frequency analysis according to claim 1 is characterized in that: For earthquake signal S EQ The way to perform wavelet transform is: According to the preset decomposition scale number n, the seismic signal is decomposed into n frequency components F(S EQ ), the decomposition formula is: In the formula, S j,k is the wavelet coefficient, which represents the local features at a given scale j and position k; S EQ (t) is the seismic signal at time t; is the wavelet function; a is the scale parameter; k is the displacement parameter; t is the time; Define the frequency component F(S EQ ) is calculated as: Where J represents the wavelet transform of the seismic signal S EQ The maximum scale for decomposition; K represents the wavelet coefficient S at each scale j j,k The number of locations; F(j) represents the frequency value corresponding to a specific scale j, Where C is a constant, a j is the scale factor associated with scale j.

3. The multi-scale and OVT denoising processing system for seismic signals based on time-frequency analysis according to claim 2, characterized in that: The high frequency components are processed in the following manner: Extract the frequency component F(S EQ ) high-frequency components, the extraction formula is: In the formula, H(F(S EQ )) represents the transfer function of the high-pass filter; F C Preset frequency value for high pass filter;. Calculate the mean and standard deviation of the high-frequency components using the following formula: Where μ is the mean value of the high-frequency component; N is the total number of signal values ​​in the high-frequency component; is the signal value of the ith high-frequency component; σ is the standard deviation of the high-frequency component; Set the dynamic artifact threshold TH ART =μ+κ·σ; where κ is a constant coefficient; The mean and standard deviation of the high-frequency components are updated based on the sliding window. The updating method is: Set the size of the sliding window and the initial mean μ init , standard deviation σ init and the corresponding dynamic artifact threshold TH ARTinit Slide the sliding window forward by 1 sample, remove the earliest sample in the sliding window, add the next sample, calculate the new mean and standard deviation, and update the dynamic artifact threshold to get the mean μ of the next sliding window init+1 , standard deviation σ init+1 and the dynamic artifact threshold TH ARTinit+1 ; Repeat the above process until all samples are processed; finally get the processed high frequency component COMP h ' igh .

4. The seismic signal multi-scale and OVT denoising processing system based on time-frequency analysis according to claim 3 is characterized in that: The calculation formula of the signal-to-noise ratio SNR of the effective signal is: Where P signal is the signal power of the effective signal, and the calculation formula is: Where Q is the number of samples of the valid signal, q is the serial index of the valid signal; P noise is the noise power, and the calculation formula is: Where P is the number of noise samples, Noise p is the remaining noise sample, obtained based on the above waveform comparison and analysis, and p is the serial index of the noise.

5. The multi-scale and OVT denoising processing system for seismic signals based on time-frequency analysis according to claim 4, characterized in that: The method of separating out the independent components includes: The signal-to-noise ratio SNR is greater than the signal-to-noise ratio threshold TH SNR The effective signal S eff Extract, center and whiten it to get the effective signal after processing After processing, the effective signal Separate and extract the independent component IC by: Set the initial weight matrix W to convert the processed valid signal Mapped to the matrix of the independent component space; Update the initial weight matrix W until it converges; the update formula is: Where W new is the new weight matrix obtained after the current iteration, W old is the old weight matrix obtained in the previous iteration; g(...) is the nonlinear activation function; g′(...) is the derivative of the nonlinear activation function; ψ{...} is the expectation operator; The formula for extracting independent components is: Where W final is the final weight matrix after iterative update.

6. The seismic signal multi-scale and OVT denoising processing system based on time-frequency analysis according to claim 5, characterized in that: The reconstruction error E is calculated re The formula is: Where M is the effective signal after processing The number of samples.

7. The multi-scale and OVT denoising processing system for seismic signals based on time-frequency analysis according to claim 5, characterized in that: The method of constructing the tide and current impact model is as follows: Based on tidal data S Tide and ocean current data S Ocea Calculate the change in signal propagation speed ΔV = a·S Tide +b·S Ocean +c; where ΔV represents the change in the propagation velocity of seismic signals under the influence of tides and ocean currents; a represents the influence of tidal changes on the propagation velocity of seismic signals, b represents the influence of ocean currents on the propagation velocity of seismic signals, and c represents the basic propagation velocity; Calculate the adjusted propagation speed V based on the signal propagation speed change ΔV adj =ΔV+V init , where V init is the initial propagation speed; Based on the adjusted propagation speed V adj For the reconstructed signal S re The wave equation is updated to obtain Where P EQ represents the pressure fluctuation of the reconstructed signal; is the Laplace operator; the pressure fluctuation P of the reconstructed signal EQ Perform Hilbert transform to obtain the complex signal of the reconstructed signal Where P EQ (t) is the pressure fluctuation at time point t; P EQ The Hilbert transform result of ; χ is an imaginary unit; the instantaneous amplitude A(t) and instantaneous phase φ(t) are extracted from the complex signal X(t) of the reconstructed signal, and the extraction formula is: A(t) = |X(t)|, φ(t) = arg(X(t)); The correction signal is calculated using the instantaneous amplitude A(t) and instantaneous phase φ(t). The calculation formula is: S Corr (t) = A(t)·e j φ(t) ; In the formula, e jφ(t) is the complex exponential form; S Corr (t) is the correction signal.

8. The seismic signal multi-scale and OVT denoising processing system based on time-frequency analysis according to claim 7, characterized in that: The training method of the tide and current impact model is as follows: Define a recursive neural network model as the basic structure of the tide and current impact model, the basic structure includes an input layer, a hidden layer, and an output layer; Introduce the parameter vector time step d, and transform the vector time step d and tidal data S Tide and ocean current data S Ocean Combined, we get the fused feature vector X d ={(S Tide ,S Ocean ) d-1 ,(S Tide ,S Ocean ) d-2 ,...,(S Tide ,S Ocean ) d-s }, where s is the number of vector time steps d; Collect historical tide data S within a fixed period of time in the past Tide and historical ocean current data S Ocean , and correspondingly construct the historical fusion feature vector Fusion of history into feature vector Expand by vector time steps to form sequence data as training input for the tide and current impact model; Initialize the network parameters of the tide and current impact model; the network parameters include the weight matrix W from the input layer to the hidden layer x , the recurrent weight matrix W of the hidden layer h , the bias vector b of the hidden layer h , the weight matrix W from the hidden layer to the output layer o and the bias vector b of the output layer o ; Define the loss function for the model affected by tides and currents In the formula, S is the length of the sequence data; Y true,d is the true value of the dth vector time step; Y pred,d is the predicted value of the dth vector time step; The sequence data is passed to the input layer, and then passes through the hidden layer and the output layer in turn, the predicted value is output, the corresponding loss function value is calculated, and the error gradient is back-propagated from the output layer to the input layer; according to the error gradient, the network parameters are updated using the optimization algorithm; the sequence data is repeatedly trained until the tide and current impact model converges or reaches the preset number of iterations, which means that the tide and current impact model training is completed.

9. The seismic signal multi-scale and OVT denoising processing system based on time-frequency analysis according to claim 8, characterized in that: The input layer is used to receive the fused feature vector X d As input; the hidden layer consists of several state vectors, and the calculation formula for each state vector is: d =f(W h ·h d-1 +W x ·X d +b h );where h d is the hidden state of the current vector at time step d, and f(...) is the activation function; The output layer calculates the output results of the tide and current impact model based on the output of the hidden layer. The calculation formula is: S Corr (t) = η(W O ·h d +b o ), where η(...) is the activation function of the output layer.

10. The seismic signal multi-scale and OVT denoising processing system based on time-frequency analysis according to claim 8, characterized in that: The method for updating the tide and current impact model is as follows: Calculate the ocean depth D at the current acquisition location Ocean The sound wave propagation speed V(D Ocean ,f EQ )=V0+ΔV(D Ocean );where V0 is the speed of sound at unit speed; ΔV(D Ocean ) is the ocean depth D at the current acquisition location Ocean The change in the speed of sound waves under water temperature, salinity and pressure is calculated as: ΔV(D Ocean )=μ1·WT+μ2·SAL+μ3·WP; where μ1 is the influence coefficient of water temperature on the sound wave velocity, μ2 is the influence coefficient of salinity on the sound wave velocity, and μ3 is the influence coefficient of water pressure on the sound wave velocity; Calculate signal propagation time Where, the signal propagation time T Prop (D Ocean ,f EQ ) represents the signal propagation time from the source to the receiver; Dist represents the distance from the source to the receiver, Where (x RX ,y RX ,z RX ) are the coordinates of the receiver; The change in the speed of the sound wave ΔV(D Ocean ) and signal propagation time T Prop (D Ocean ,f EQ ) is returned to the tide and current influence model to correct the signal S Corr (t) is updated to obtain the enhanced signal S Reinf (T(D Ocean ,f EQ )).

Images