Surge Correction Method, System and Storage Medium Based on Variational Mode Decomposition

Abstract

The present invention discloses a surge correction method, system and storage medium based on variational mode decomposition, which relates to the technical field of surge correction. The specific steps are as follows: obtaining shallow stratum profile data and picking up the arrival time curve of the seabed reflection wave; extracting the seabed topography trend term from the arrival time curve of the seabed reflection wave through variational mode decomposition; calculating the surge correction amount on each trace according to the arrival time curve of the seabed reflection wave and the seabed topography trend term, and performing surge correction through the surge correction amount. By replacing the traditional filtering method with variational mode decomposition, the present invention adaptively extracts the topography trend from the seabed reflection isochrone curve, effectively solves the problem of isochrone jitter of the shallow stratum profile caused by surges, has strong anti-noise ability, high correction accuracy, does not require frequent manual estimation of surge parameters, and is suitable for batch processing of large-scale data.

Description

Technical Field

[0001] The present invention relates to the technical field of surge correction, and particularly to a surge correction method, system and storage medium based on variational mode decomposition. Background Art

[0002] A shallow seismic profiler is a geophysical technical means for continuously detecting information such as shallow subsurface stratigraphic structures, tectonics, and disaster geology underwater based on acoustic principles. With characteristics such as high vertical resolution and efficient operation in shallow seismic profiling, it is widely used in marine engineering geological surveys, late Quaternary stratigraphy research, buried object detection, port construction and daily maintenance. When conducting shallow seismic profiling, sea waves or surges can affect the attitude of the ship or the carrier of the shallow seismic profiler, thereby affecting the quality of the shallow seismic profile. When the survey ship is sailing, its draft depth remains basically unchanged, but its elevation fluctuates due to sea waves or surges, resulting in a change in the actual propagation path of the seismic wave excited by the transducer compared to the theoretical path. This change in propagation distance causes a regular dislocation of the isochrones of the seabed and underlying strata in the shallow seismic profile. A typical feature is that the isochrone of the seabed shows an up-and-down jitter phenomenon, thus significantly reducing the signal-to-noise ratio and resolution of the profile. Therefore, how to effectively reduce or even eliminate this influence is the key to obtaining a high-quality shallow seismic profile, and this process is also called surge correction.

[0003] In the prior art, methods for surge correction include: surge correction based on the cross-correlation method of model traces; picking up the seabed and then filtering, and taking the difference in reflection time before and after filtering as the correction amount to eliminate the influence of surges; using external constraint information such as the seabed depth derived from SBP and MBES for surge correction, etc. Surge correction with the aid of external information (water depth data) strongly depends on the availability of measured water depth data in the survey area and the accuracy of the water depth data, and its application is restricted. The method based on model traces requires aligning the signals in the adjacent area of the model trace, and the model trace needs to be frequently updated during calculation, resulting in low operation efficiency. Currently, the method of picking up the seabed and then filtering has been widely used, but the parameter setting of the filter will affect the effect of surge correction, and it is necessary to manually estimate the period of the surge to match the filter parameters. The periods of surges under different sea conditions are not consistent. In order to adaptively process a large number of shallow profile data, a surge correction method that can be adaptive and less affected by parameters is needed. Summary of the Invention

[0004] The purpose of the present invention is to provide a surge correction method, system and storage medium based on variational mode decomposition to solve the problems raised in the background art.

[0005] To achieve the above purpose, the present invention provides the following solution: A surge correction method based on variational mode decomposition, and the specific steps include the following:

[0006] Obtain shallow stratum profile data and pick up the arrival time curve of the seabed reflection wave;

[0007] Extract the seabed terrain trend term from the arrival time curve of the seabed reflection wave through variational mode decomposition;

[0008] Calculate the swell correction amount on each channel according to the arrival time curve of the seabed reflection wave and the seabed terrain trend term, and perform swell correction through the swell correction amount.

[0009] Preferably, read the shallow stratum profile data and pick up the arrival time curve of the seabed reflection wave according to the quality of the profile by using the first arrival picking method.

[0010] Preferably, the specific process of extracting the seabed terrain trend term through the variational mode decomposition is as follows:

[0011] Assume that the arrival time curve U1(t) of the seabed reflection wave is composed of frequency modulation - amplitude modulation signals, and the expression is: Among them, the phase is a non - decreasing function, and A(t) is the envelope;

[0012] Use u k (t) to represent each obtained mode function, and perform Hilbert transform on each mode function to obtain its analytic signal: Among them, δ(t) is the impulse function; is the singular function, which is used to construct the analytic signal and represent its frequency components;

[0013] Estimate the center frequency of each mode analytic signal, modulate each mode spectrum to the corresponding base frequency band, and obtain the corresponding constrained variational problem:

[0014] Solve the constrained variational problem and update to find the optimal solution of the constrained variational problem.

[0015] Preferably, the expression of the constrained variational problem is:

[0016]

[0017] Among them, {u k} = {u1,..., u K} represents k modes, and each mode is a modal; {ω k} = {ω1,..., ω K} represents the center frequency corresponding to each mode, is the center frequency of each mode analytic signal;

[0018] Introduce L2 regularization and Lagrange multipliers to solve the constrained variational problem and update to find the optimal solution of the constrained variational problem; update to find the optimal solution of the constrained variational problem. Successively obtain u k (t), and the corresponding center frequency is ω k , from high frequency to low frequency, select the mode with the lowest center frequency as the seabed terrain trend, where k is the number of modes selected during decomposition.

[0019] Preferably, updating to find the optimal solution of the constrained variational problem includes mode update, center frequency update, and dual ascent update.

[0020] Preferably, subtract the seabed reflection wave arrival time curve from the seabed terrain trend term to obtain the swell correction amount for each trace.

[0021] On the other hand, provide a swell correction system based on variational mode decomposition, including a data acquisition module, a data analysis module, and a correction module; wherein,

[0022] The data acquisition module is used to acquire shallow stratigraphic profile data and pick up the seabed reflection wave arrival time curve;

[0023] The data analysis module is used to extract the seabed terrain trend term from the seabed reflection wave arrival time curve through variational mode decomposition;

[0024] The correction module is used to calculate the swell correction amount for each trace according to the seabed reflection wave arrival time curve and the seabed terrain trend term, and perform swell correction through the swell correction amount.

[0025] Finally, provide a computer storage medium, on which a computer program is stored, and when the computer program is executed by a processor, the steps of the swell correction method based on variational mode decomposition are implemented.

[0026] According to the specific embodiments provided by the present invention, the following technical effects are disclosed by the present invention:

[0027] (1) The present invention can perform swell correction adaptively, without the need to pre-estimate the influence period of the swell, with little influence from parameter selection, and is suitable for processing a large number of shallow profile data. After processing, the seabed reflection isochrones on the shallow stratigraphic profile are flattened, and the underlying sedimentary layers are also clearly visible, improving the profile quality.

[0028] (2) By applying the VMD method instead of the traditional filtering method, the present invention adaptively extracts the terrain trend from the seabed reflection isochrone curve, effectively solves the problem of isochrone jitter of the shallow stratigraphic profile caused by swell, has strong anti-noise ability, high correction accuracy, and does not require frequent manual estimation of swell parameters, and is suitable for batch processing of large-scale data. Description of the Drawings

[0029] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the accompanying drawings required in the embodiments. Obviously, the accompanying drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other accompanying drawings can be obtained based on these drawings.

[0030] Figure 1 It is the flowchart of the method of the present invention;

[0031] Figure 2 It is the arrival time curve graph of the picked seabed reflection wave of the present invention;

[0032] Figure 3 It is the ultra-high frequency interference graph obtained by VMD decomposition of the present invention;

[0033] Figure 4 It is the schematic diagram of IMF2 obtained by VMD decomposition of the present invention;

[0034] Figure 5 It is the schematic diagram of the seabed terrain trend term obtained by VMD decomposition of the present invention;

[0035] Figure 6 It is the comparison graph of the shallow profile images before and after flat terrain correction of the present invention;

[0036] Figure 7 It is the comparison graph of the shallow profile images before and after rugged terrain correction of the present invention;

[0037] Figure 8 It is the seabed arrival time curve graph after superimposing the wave height of 15 and the wavelength of 60 of the present invention;

[0038] Figure 9 It is the seabed arrival time curve graph after superimposing the wave height of 15 and the wavelength of 80 of the present invention. Detailed implementation manners

[0039] The following will clearly and completely describe the technical solutions in the embodiments of the present invention with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, rather than all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts belong to the scope of protection of the present invention.

[0040] The object of the present invention is to provide a surge correction method based on variational mode decomposition, as Figure 1 shown, and the specific steps are as follows:

[0041] S1. Obtain shallow stratum profile data and pick up the arrival time curve of the seabed reflection wave;

[0042] S2. Extract the seabed terrain trend term from the seabed reflection wave arrival time curve through variational mode decomposition;

[0043] S3. Calculate the swell correction amount for each trace according to the seabed reflection wave arrival time curve and the seabed terrain trend term, and perform swell correction through the swell correction amount.

[0044] The seabed reflection wave arrival time is the two-way travel time when the seismic wave excited by the sub-bottom profiler is reflected by the seabed and the echo is first recorded by the sub-bottom profiler; on the seismic shot gather, the seabed reflection wave arrival times corresponding to each shot are combined into a seabed reflection wave arrival time curve.

[0045] Further, step S1 specifically includes:

[0046] Assume that there are N traces in the shallow stratum profile data, and each trace has m sampling points. The seismic recording amplitude of the i-th trace is S i (t), the obtained seabed reflection wave arrival time is Tbi, and the entire profile is combined into an arrival time curve U1 = [Tb1, Tb2,..., Tbn].

[0047] Seabed picking in shallow stratum profile data processing is similar to first arrival picking of seismic data. Seismic first arrival picking methods can be divided into traditional methods and machine learning methods in research; traditional methods mainly include the energy ratio method, the AIC method, the fractal dimension method, and other methods; machine learning methods mainly include clustering algorithm-based, classification algorithm-based, and image classification or segmentation methods. Usually, shallow profile observation data is presented in the form of common receiver gather. Due to the strong seabed reflection signal, in most cases, the energy ratio method can be used. The energy ratio method has a simple principle and wide application, but it is sensitive to noise. It is particularly important to select a suitable seabed picking method according to the actual situation of the data. The seabed arrival time curve obtained by using the traditional energy ratio method for actual shallow profile data in a certain area is as Figure 2 shown.

[0048] Further, step S2 specifically includes:

[0049] Perform variational mode decomposition (VMD) on the seabed reflection wave arrival time curve U1 obtained in the above step S1 to obtain each order mode (IMF1, 2,..., k) of the corresponding curve, where k is the number of modes selected during decomposition.

[0050] Variational Mode Decomposition (VMD) is a novel signal processing method. Its principle is to decompose a signal into a set of Intrinsic Mode Functions (IFMs) with specific sparse properties. The IMF is defined as an amplitude-modulated and frequency-modulated (AM-FM) signal with narrowband characteristics. This method uses Wiener filtering to directly update the modes in the Fourier domain, where the filter is tuned to the current center frequency, which is used as the centroid of the power spectrum of the mode. VMD decomposition has strong noise robustness and a solid mathematical theory foundation, is applicable to non-stationary and non-linear signal analysis, and has been widely used in various fields.

[0051] The shape of the arrival time curve of the seabed reflection wave in a shallow seismic profile is affected by various factors such as the seabed morphology and swell. Among them, the influence caused by the swell is a relatively stationary high-frequency signal, while the seabed morphology can be regarded as a non-linear low-frequency signal. The seabed reflection isochron in the profile can be considered as the result of the superposition of the swell signal and the seabed morphology trend term. The VMD method is an adaptive and quasi-orthogonal signal decomposition method, which can extract the trend term from the signal without prior information. And if the center frequency of the power spectrum of the trend term plays a dominant role in the power spectrum of the original signal, the first-order IMF component extracted in the VMD method is the trend term of the signal and is independent of the number of modes. Moreover, the seabed morphology is the main component of the arrival time curve compared with the undulation of the sea surface. Therefore, the VMD method can be used to adaptively extract the seabed morphology as the trend term from the picked arrival time curve of the seabed reflection wave. Then, by subtracting the arrival time corresponding to the seabed morphology from the picked arrival time of the reflection wave, the swell correction amount can be obtained.

[0052] The general idea of performing variational mode decomposition is as follows: for the arrival time curve U1(t) of the seabed reflection wave, find its corresponding k IMFs, so that the sum of the estimated bandwidths of each IMF around its center frequency is minimized, and ensure that the sum of all IMFs is the input arrival time curve U1(t).

[0053] Assume that the arrival time curve U1(t) of the seabed reflection wave is composed of frequency-modulated and amplitude-modulated signals, and its expression is: where the phase is a non-decreasing function, and A(t) is the envelope;

[0054] Use u k (t) to represent the obtained modal functions, and perform Hilbert transform on each modal function to obtain its analytic signal: where δ(t) is the impulse function;

[0055] Estimate the center frequency of each modal analytic signal Modulate each modal spectrum to the corresponding base frequency band:

[0056] The corresponding constrained variational problem is obtained as follows:

[0057]

[0058] where {u k} = {u1, …, u K} represents k modes, and each mode is a modality; {ω k} = {ω1, …, ω K} represents the center frequency corresponding to each mode, is the center frequency of the analytic signal of each modality;

[0059] The bandwidth of each mode is determined by the square of the H1 norm of its baseband-offset Hilbert complementary analytic signal (only with positive frequencies). Quadratic penalty and Lagrange multipliers are introduced to solve the problem in Equation (1). The augmented Lagrangian is set using the following equation:

[0060]

[0061] where α is the balance parameter for the data fidelity constraint, λ is the Lagrange multiplier, and f represents the seabed reflection arrival-time curve u(t).

[0062] The alternating direction method of multipliers (ADMM) is used to solve the variational problem in Equation (2), and different decomposed modes and center frequencies are generated during each shift operation. Each modality obtained through the solution in the frequency domain can be expressed as:

[0063]

[0064] where is the instantaneous frequency, and ω represents the center frequency of the modality.

[0065] The optimal solution of the variational problem is found by updating , where the superscript n+1 represents the result of the (n + 1)-th iteration, and the update process is as follows:

[0066] 1. Mode update: The update of mode is as shown in Equation (4). Embedding Wiener filtering, the mode is directly updated in the Fourier domain through a filter tuned to the current center frequency :

[0067]

[0068] (2 - 5)

[0069] 2. Center frequency update: The center frequency is updated to the centroid of the power spectrum of the corresponding modality, as shown in Equation (5):

[0070]

[0071] 3. Dual ascent update: For all ω ≥ 0, update the Lagrange multipliers by Equation (6). Dual ascent enforces exact signal reconstruction until

[0072]

[0073] The result obtained after picking up the arrival time curve from the actual data and performing VMD decomposition is as follows Figure 2-4 shown. In this embodiment, k = 3 is taken during decomposition (the present invention is not sensitive to the selection of parameters. When the picked surge noise is large, it is recommended to select a K value greater than 10. When the picked surge noise is small, the K value can be selected smaller, such as 3, to improve the operation speed). Among them Figure 2 is the picked arrival time curve of the seabed reflection wave. There are up and down jitters with inconsistent periods and magnitudes on this curve, which are caused by surges; Figure 3 is IMF1, which is the ultra-high frequency interference obtained by VMD decomposition and may be caused by the error in arrival time extraction and recording error; Figure 4 is IMF2, which is the influence of the surge with a certain periodicity obtained by VMD decomposition on the arrival time; Figure 5 is IMF3. Among them, IMF3 is the ultra-low frequency trend term in the arrival time curve. It is considered that this IMF is the trend term of the seabed morphology, and the seabed terrain trend term U2 is obtained, where U2 = [Ts1, Ts2,..., Tsn].

[0074] Furthermore, step S3 specifically includes:

[0075] Based on the picked arrival time curve U1 of the seabed reflection wave and the seabed terrain trend term U2 obtained in the previous step S2, subtract U1 from U2 to obtain the time shift amount ΔT = [ΔT1, ΔT2,..., ΔTn] of each channel due to the arrival time of the surge. When ΔTi > 0, it indicates that the surge causes the position of the profiler in the current observation channel of the profile to rise relative to the reference plane, resulting in an increase in the two-way travel time of the seabed reflection wave. Then, the data of this channel needs to be time-shifted upward as a whole. On the contrary, when ΔTi < 0, the data of this channel is time-shifted downward.

[0076] The process of time shift is as follows:

[0077] tci = tbi - Δti;

[0078] Δpoint = ceil(tci / Tslice);

[0079] The time before correction for each trace is tbi, and the time after correction is tci. Δpoint represents the number of sampling points shifted up or down. When Δpoint > 0, it means removing the first Δpoint sampling points of the data for that trace and appending Δpoint zeros at the end. Similarly, when Δpoint < 0, it means appending Δpoint zeros at the beginning of the data for that trace and removing the last Δpoint sampling points of the data. Here, ceil represents rounding down, and Tslice represents the sampling interval of the shallow profile record.

[0080] The effect diagrams before and after using VMD decomposition correction for the actual data of a certain place are as Figure 6-7 shown Figure 6 The comparison of the shallow profile images before and after flat terrain correction is Figure 7 The comparison of the shallow profile images before and after rugged terrain correction. It can be seen that after correction using the method of the present invention, whether it is flat terrain or rugged terrain, the seabed no longer jitters, the strata are clearly visible, and the continuity is enhanced.

[0081] After generating the seabed terrain trend term through forward modeling simulation and superimposing the surf signal to simulate noise on it, use the VMD method mentioned in the present invention and other filtering methods to process the arrival time curve of the seabed reflection wave after superposition respectively. Calculate the correlation coefficient between the terrain trend obtained after processing and the seabed terrain trend term generated by the original forward modeling to measure the effect after processing. The closer the correlation coefficient is to 1, the closer the terrain trend term obtained after processing is to the preset terrain trend, that is, the better the effect. The results obtained are as Figure 8-9 shown. It is not difficult to see that the VMD method has better accuracy and anti-noise performance compared with the traditional filter method, and its results are stable and less affected by parameters.

[0082] On the other hand, an embodiment of the present invention provides a surf correction system based on variational mode decomposition, including a data acquisition module, a data analysis module, and a correction module; wherein,

[0083] The data acquisition module is used to acquire shallow stratigraphic profile data and pick up the arrival time curve of the seabed reflection wave;

[0084] The data analysis module is used to extract the seabed terrain trend term from the arrival time curve of the seabed reflection wave through variational mode decomposition;

[0085] The correction module is used to calculate the surf correction amount for each trace according to the arrival time curve of the seabed reflection wave and the seabed terrain trend term, and perform surf correction through the surf correction amount.

[0086] Finally, a computer storage medium is provided. A computer program is stored on the computer storage medium. When the computer program is executed by a processor, it realizes the steps of a surf correction method based on variational mode decomposition.

[0087] In this article, specific examples are used to elaborate on the principles and implementation manners of the present invention. The description of the above embodiments is only used to help understand the method and its core idea of the present invention. At the same time, for those of ordinary skill in the art, according to the idea of the present invention, there will be changes in the specific implementation manners and application scopes. In summary, the content of this specification should not be construed as a limitation to the present invention.

Claims

1. A surge correction method based on variational mode decomposition, characterized in that, The specific steps are as follows: Obtain shallow stratum profile data and pick up the arrival time curve of the seabed reflection wave; Extract the seabed topography trend term from the arrival time curve of the seabed reflection wave through variational mode decomposition; Calculate the swell correction amount on each trace according to the arrival time curve of the seabed reflection wave and the seabed topography trend term, and perform swell correction through the swell correction amount; Read the shallow stratum profile data and pick up the arrival time curve of the seabed reflection wave according to the quality of the profile by using the first arrival picking method; Extract the seabed topography trend term from the arrival time curve of the seabed reflection wave by using the variational mode decomposition, and the specific process is as follows: Assume that the arrival time curve U1(t) of the seabed reflection wave is composed of frequency modulation - amplitude modulation signals, and the expression is: where the phase is a non-decreasing function and A(t) is the envelope; The obtained modal functions are represented by u k (t). Hilbert transform is performed on each modal function to obtain its analytic signal: where, δ(t) is the impulse function; is a singular function used to construct the analytic signal and represent its frequency components; Estimate the center frequency of each modal analytic signal, modulate each modal spectrum to the corresponding base frequency band, and obtain the corresponding constrained variational problem: where, {u k} = {u1, …, u k} represents k patterns, and each pattern is a modality; {ω k} = {ω1, …, ω k} represents the central frequency corresponding to each pattern, being the central frequency of the analytic signal of each modality; Introduce L2 regularization and Lagrange multipliers to solve the constrained variational problem and update to find the optimal solution of the constrained variational problem; by updating to find the optimal solution of the constrained variational problem, u k (t) is obtained in sequence, and the corresponding center frequency is ω k , from high frequency to low frequency, select the mode with the lowest center frequency as the seabed topography trend, where k is the number of modes selected during decomposition.

2. The surge correction method based on variational mode decomposition according to claim 1, wherein, Subtract the seabed topography trend term from the arrival time curve of the seabed reflection wave to obtain the swell correction amount on each trace.

3. A surge correction system based on variational mode decomposition, characterized in that, It is used to implement a swell correction method based on variational mode decomposition as described in claim 1, including a data picking module, a data analysis module, and a correction module; wherein, The data picking module is used to obtain shallow stratum profile data and pick up the arrival time curve of the seabed reflection wave; The data analysis module is used to extract the seabed topography trend term from the arrival time curve of the seabed reflection wave through variational mode decomposition; The correction module is used to calculate the swell correction amount on each trace according to the arrival time curve of the seabed reflection wave and the seabed topography trend term, and perform swell correction through the swell correction amount.

4. A computer storage medium, characterized in that, A computer program is stored on the computer storage medium, and when the computer program is executed by a processor, it implements the steps of a swell correction method based on variational mode decomposition as described in any one of claims 1 - 2.

Images