Surge correction method and system based on variational mode decomposition, and storage medium

Abstract

The invention discloses a surge correction method and system based on variational mode decomposition, and a storage medium, and relates to the technical field of surge correction. The method comprises the following specific steps: acquiring shallow stratum profile data, and picking up a seabed reflected wave arrival time curve; extracting submarine topography trend terms from the submarine reflected wave arrival time curve through variational mode decomposition; and according to the submarine reflected wave arrival time curve and the submarine topography trend term, calculating the surge correction quantity on each channel, and carrying out surge correction through the surge correction quantity. According to the method, a traditional filtering method is replaced by variational mode decomposition, the terrain trend is adaptively extracted from the seabed reflection event curve, the problem of shallow stratum profile event jitter caused by the surge is effectively solved, the anti-noise capability is high, the correction precision is high, frequent manual estimation of surge parameters is not needed, and the method is suitable for large-scale popularization and application. The method 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 in particular to a surge correction method, system and storage medium based on variational mode decomposition. Background Art

[0002] The shallow layer profiler is a geophysical technology based on acoustic principles to continuously detect underwater shallow stratum structure, structure, disaster geology and other information. Shallow layer profile detection has the characteristics of high vertical resolution and efficient operation, and is widely used in marine engineering geological surveys, late Quaternary stratum research, buried object detection, port construction and daily maintenance. When measuring shallow layer profiles, waves or swells will affect the posture of the ship or the carrier of the shallow layer profiler, thereby affecting the quality of the shallow layer profile. The draft of the survey ship remains basically unchanged when it is sailing, but its elevation will fluctuate due to waves or swells, thereby changing the actual propagation path and theoretical path of the seismic waves excited by the transducer. This change in propagation distance will cause the phase axis of the seabed and the underlying strata in the shallow layer profile to shift regularly. The typical feature is that the phase axis of the seabed shows up and down shaking, which significantly reduces 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 high-quality shallow layer profiles. This process is also called swell correction.

[0003] The methods for surge correction in the prior art include: surge correction based on the model channel cross-correlation method; filtering after picking up the seabed, using the difference in reflection time before and after filtering as the correction amount to eliminate the influence of surge; using external constraint information such as SBP, MBES derived seabed depth for surge correction, etc. Surge correction with the help of external information (water depth data) strongly depends on the presence of measured water depth data and the accuracy of water depth data in the survey area, and its application is restricted. The method based on the model channel requires the alignment of the signals in the adjacent area of ​​the model channel, and the model channel needs to be updated frequently during calculation, resulting in low computational efficiency. At present, the method of filtering after picking up the seabed has been widely used, but the parameter setting of the filter will affect the effect of surge correction, and it is necessary to realize artificial estimation of the surge period to match the filter parameters. The period of surges is not consistent under different sea conditions. In order to adaptively process a large number of shallow profile data, a surge correction method that can be adaptive and is 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 modal decomposition to solve the problems raised in the background technology.

[0005] To achieve the above object, the present invention provides the following solution: a surge correction method based on variational mode decomposition, the specific steps of which include the following:

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

[0007] Extracting a seafloor topography trend term from the seafloor reflection wave arrival time curve by variational mode decomposition;

[0008] The surge correction amount on each channel is calculated according to the arrival time curve of the seabed reflection wave and the seabed topography trend item, and the surge correction is performed by using the surge correction amount.

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

[0010] Preferably, the specific process of extracting the seabed topography trend term by the variational mode decomposition is:

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

[0012] Use u to obtain each modal function k (t) represents that the Hilbert transform is performed on each modal function to obtain its analytical signal: Among them, δ(t) is the pulse function; is a singular function used to construct an analytical signal and represent its frequency components;

[0013] Estimate the center frequency of each modal analytical signal, modulate each modal spectrum to the corresponding baseband, and obtain the corresponding constrained variational problem:

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

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

[0016]

[0017] Among them, {u k}={u1,…,u K} represents k modes, each of which is a mode; {ω k}={ω1,…,ω K} represents the center frequency corresponding to each mode, is the center frequency of each modal analytical 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; find the optimal solution of the constrained variational problem by updating. Get u in turn k (t), the corresponding center frequency is ω k , from high frequency to low frequency, the mode with the lowest central frequency is selected as the seabed topography trend, where k is the number of modes selected during decomposition.

[0019] Preferably, the optimal solution of the constrained variational problem is found through updates, including mode updates, center frequency updates, and dual ascent updates.

[0020] Preferably, the seabed reflection wave arrival time curve is subtracted from the seabed topography trend term to obtain the surge correction amount on each channel.

[0021] On the other hand, a surge correction system based on variational mode decomposition is provided, comprising a data pickup module, a data analysis module, and a correction module; wherein:

[0022] The data picking module is used to obtain shallow stratum profile data and pick up the arrival time curve of the seabed reflection wave;

[0023] The data analysis module is used to extract the seabed topography trend item by variational mode decomposition of the seabed reflection wave arrival time curve;

[0024] The correction module is used to calculate the surge correction amount on each channel according to the seabed reflection wave arrival time curve and the seabed topography trend item, and perform surge correction using the surge correction amount.

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

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

[0027] (1) The present invention can adaptively perform surge correction without pre-estimating the surge impact period. The influence of parameter selection is small, and it is suitable for processing a large number of shallow profile data. After processing, the seabed reflection event axis on the shallow stratum profile is flattened, and the underlying sedimentary layer is also clearly visible, and the profile quality is improved.

[0028] (2) The present invention uses the VMD method instead of the traditional filtering method to adaptively extract the terrain trend from the seabed reflection event curve, effectively solving the problem of shallow layer profile event jitter caused by surges. It has strong anti-noise ability and high correction accuracy, and does not require frequent manual estimation of surge parameters. It is suitable for batch processing of large-scale data. BRIEF DESCRIPTION OF THE DRAWINGS

[0029] 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 will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0030] Figure 1 is a flow chart of the method of the present invention;

[0031] Figure 2 The picked-up seabed reflection wave arrival time curve diagram of the present invention;

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

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

[0034] Figure 5 A schematic diagram of the seafloor topography trend item obtained by VMD decomposition of the present invention;

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

[0036] Figure 7 A comparison diagram of the shallow profile image before and after the correction of the rough terrain of the present invention;

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

[0038] Fig. 9 This is a seabed arrival time curve diagram after superimposing a wave height of 15 and a wavelength of 80 according to the present invention. DETAILED DESCRIPTION

[0039] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only 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.

[0040] The purpose of the present invention is to provide a surge correction method based on variational mode decomposition, such as Figure 1 As shown, the specific steps include the following:

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

[0042] S2, extracting the seabed topography trend term from the seabed reflection wave arrival time curve by variational mode decomposition;

[0043] S3. Calculate the surge correction on each channel based on the arrival time curve of the seabed reflection wave and the seabed topography trend item, and perform surge correction using the surge correction.

[0044] The arrival time of the seabed reflection wave is the two-way travel time of the echo first recorded by the profiler after the seismic wave excited by the profiler is reflected by the seabed; in the seismic shot collection, the arrival time of the seabed reflection wave corresponding to each shot is combined into a seabed reflection wave arrival time curve.

[0045] Furthermore, step S1 specifically includes:

[0046] Assume that there are N shallow stratum profile data, each with m sampling points, and the amplitude of the seismic record of the ith track is S i (t), the arrival time of the picked-up seabed reflection wave is Tbi, and the entire profile is combined into an arrival curve U1 = [Tb1, Tb2, ..., Tbn].

[0047] Seabed picking in shallow 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 energy ratio method, AIC method, fractal dimension method and other methods; machine learning methods mainly include clustering algorithm-based, classification algorithm-based and image classification or segmentation methods. Normally, shallow profile observation data are presented in the form of common receiving point gathers. Due to the strong seabed reflection signal, the energy ratio method can be used in most cases. The energy ratio method has a simple principle and is widely used, but it is sensitive to noise. It is particularly important to choose a suitable seabed picking method based on the actual situation of the data. The seabed arrival time curve obtained by using the traditional energy ratio method to pick the seabed for actual shallow profile data in a certain area is shown below. Figure 2 shown.

[0048] Furthermore, step S2 specifically includes:

[0049] The seabed reflection wave arrival curve U1 obtained in the above step S1 is subjected to variational mode decomposition (VMD) to obtain the modes of each order (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 the signal into a set of intrinsic mode functions (IFMs) with specific sparse properties. IMFs are defined as amplitude modulated frequency modulated (AM-FM) signals with narrowband characteristics. This method uses Wiener filtering to directly update the mode in the Fourier domain, where the filter is tuned to the current center frequency as the center of gravity of the modal power spectrum. VMD decomposition has strong noise robustness and a solid mathematical theoretical foundation. It is suitable for non-stationary and nonlinear signal analysis and has been widely used in various fields.

[0051] The shape of the arrival curve of the seabed reflection wave in the shallow layer profile is affected by many factors such as the seabed shape and swell. The impact caused by the swell is a relatively stable high-frequency signal, while the seabed shape can be regarded as a nonlinear low-frequency signal. The seabed reflection event axis in the profile can be considered as the result of the superposition of the swell signal and the trend term of the seabed shape. The VMD method is an adaptive, quasi-orthogonal signal decomposition method that can extract the trend term from the signal without prior information. If the power spectrum center frequency of the trend term plays a dominant role in the power spectrum of the original signal, the first-order IMF component extracted by the VMD method is the trend term of the signal, and it is independent of the value of the number of modes. The seabed shape is the main component of the arrival curve compared to the undulation of the sea surface. Therefore, the VMD method can be used to adaptively extract the seabed shape as a trend term from the picked seabed reflection wave arrival curve. Then the arrival time corresponding to the seabed shape is subtracted from the picked reflection wave arrival time to obtain the swell correction amount.

[0052] The general idea of ​​variational mode decomposition is to find the corresponding k IMFs for the arrival curve U1(t) of the seabed reflection wave, so that the sum of the estimated bandwidth of each IMF around its center frequency is minimized, and the sum of each IMF is guaranteed to be the input arrival curve U1(t).

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

[0054] Use u to obtain each modal function k (t) represents that the Hilbert transform is performed on each modal function to obtain its analytical signal: Among them, δ(t) is the pulse function;

[0055] Estimate the center frequency of each modal analytical signal Modulate each modal spectrum to the corresponding baseband:

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

[0057]

[0058] Among them, {u k}={u1,…,u K} represents k modes, each of which is a mode; {ω k}={ω1,…,ω K} represents the center frequency corresponding to each mode, is the center frequency of each modal analytical signal;

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

[0060]

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

[0062] The variational problem in equation (2) is solved using the alternating direction method of multipliers (ADMM) method to generate different decomposed modes and center frequencies during each shift operation. Through the solution in the frequency domain, each mode obtained can be expressed as:

[0063]

[0064] in, is the instantaneous frequency, and ω represents the center frequency of the mode.

[0065] By updating Find the optimal solution to the variational problem, where the superscript n+1 Represents the result of the n+1th iteration. The updating process is as follows:

[0066] 1. Mode update: Mode The update of is shown in equation (4). Embedded Wiener filter, by tuning to the current center frequency The filter updates the pattern directly in the Fourier domain:

[0067]

[0068] (2-5)

[0069] 2. Center frequency update: Center frequency Updated to the center of gravity of the corresponding modal power spectrum, as shown in formula (5):

[0070]

[0071] 3. Dual ascending update: For all ω ≥ 0, update the Lagrange multiplier by equation (6): The duality rises, forcing exact signal reconstruction, until

[0072]

[0073] The results obtained after picking up the arrival time curve from the actual data and performing VMD decomposition are as follows Figure 2-4 As shown, in this embodiment, k is 3 when decomposing (the present invention is not sensitive to the selection of parameters. When the picked-up surge noise is large, it is recommended that the K value be greater than 10. When the picked-up surge noise is small, the K value can be smaller, such as 3, to increase the calculation speed). Figure 2 This is the arrival time curve of the picked-up seabed reflection wave. The up and down jitters with inconsistent periods and sizes on this curve are caused by swells. Figure 3 IMF1 is the ultra-high frequency interference obtained by VMD decomposition, which may be caused by the error of extraction and recording; Figure 4 IMF2 is the impact of the surge with a certain periodicity on the arrival time obtained by VMD decomposition; Figure 5 IMF3 is the ultra-low frequency trend item in the arrival curve. The IMF is considered to be the trend item of the seabed morphology, and the seabed topography trend item U2 is obtained, U2 = [Ts1, Ts2, ..., Tsn].

[0074] Furthermore, step S3 specifically includes:

[0075] According to the seabed reflection wave arrival time curve U1 and the seabed topography trend item U2 obtained in the above step S2, the time shift of each channel due to the arrival of the surge wave is obtained by subtracting U1 and U2, and obtaining the time shift of each channel due to the arrival of the surge wave ΔT=[ΔT1, ΔT2,……, ΔTn]. When ΔTi is greater than zero, it indicates that the surge wave causes the position of the profiler to rise relative to the reference surface when recording the current observation channel in the profile, resulting in an increase in the two-way travel time of the seabed reflection wave, and it is necessary to shift the data of this channel upward as a whole. On the contrary, when ΔTi is less than zero, the data of this channel is shifted downward.

[0076] The time shift process is:

[0077] tci=tbi-Δti;

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

[0079] The time before correction of each channel is tbi, and the time after correction is tci. Δpoint indicates the number of sampling points to be moved up or down. When Δpoint>0, it means to remove the data of Δpoint sampling points at the beginning of the data of this channel, and add Δpoint zeros at the end. Similarly, when Δpoint<0, it means to add Δpoint sampling points of zeros at the beginning of the data of this channel, and remove Δpoint sampling points of data at the end. Among them, ceil indicates rounding down, and Tslice indicates the sampling interval of shallow section record.

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

[0081] After the seabed topography trend term is generated by forward modeling, the surge signal simulation noise is superimposed on it, and the VMD method and other filtering methods mentioned in the present invention are used to process the superimposed seabed reflection wave arrival curve respectively, and the correlation coefficient between the terrain trend obtained after processing and the seabed topography trend term generated by the original forward modeling is calculated 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 are as follows: Figure 8-9 As shown, it is not difficult to see that the VMD method has better accuracy and noise resistance than the traditional filter method, and its results are stable with less influence of parameters.

[0082] Another aspect of the present invention provides a surge correction system based on variational mode decomposition, including a data pickup module, a data analysis module, and a correction module; wherein:

[0083] Data picking module, used to obtain shallow stratum profile data and pick up the arrival time curve of seabed reflection wave;

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

[0085] The correction module is used to calculate the surge correction on each channel according to the arrival time curve of the seabed reflection wave and the trend item of the seabed topography, and perform surge correction through the surge correction.

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

[0087] The principles and implementation methods of the present invention are described in this article using specific examples. The description of the above embodiments is only used to help understand the method and core idea of ​​the present invention. At the same time, for those skilled in the art, according to the idea of ​​the present invention, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as limiting the present invention.

Claims

1. A surge correction method based on variational mode decomposition, characterized in that: The specific steps include the following: Obtain shallow stratum profile data and pick up the arrival time curve of seabed reflection waves; Extracting a seafloor topography trend term from the seafloor reflection wave arrival time curve by variational mode decomposition; The surge correction amount on each channel is calculated according to the arrival time curve of the seabed reflection wave and the seabed topography trend item, and the surge correction is performed by using the surge correction amount.

2. A surge correction method based on variational mode decomposition according to claim 1, characterized in that: The shallow stratum profile data is read, and a first arrival picking method is selected according to the quality of the profile to pick up the seabed reflection wave arrival time curve.

3. The surge correction method based on variational mode decomposition according to claim 1, characterized in that: The variational mode decomposition is used to extract the seabed topography trend term from the seabed reflection wave arrival time curve, and the specific process is as follows: Assuming that the arrival curve U1(t) of the seabed reflection wave is composed of a frequency-amplitude modulation signal, the expression is: Among them, the phase is a non-decreasing function, A(t) is the envelope; Use u to obtain each modal function k (t) represents that the Hilbert transform is performed on each modal function to obtain its analytical signal: Among them, δ(t) is the pulse function; is a singular function used to construct an analytical signal and represent its frequency components; Estimate the center frequency of each modal analytical signal, modulate each modal spectrum to the corresponding baseband, and obtain the corresponding constrained variational problem: Among them, {u k }={u1,…,u K } represents k modes, each of which is a mode; {ω k }={ω1,…,ω K } represents the center frequency corresponding to each mode, is the center frequency of each modal analytical signal; Introduce L2 regularization and Lagrange multipliers to solve the constrained variational problem and update to find the optimal solution of the constrained variational problem; find the optimal solution of the constrained variational problem by updating. Get u in turn k (t), the corresponding center frequency is ω k , from high frequency to low frequency, the mode with the lowest central frequency is selected as the seabed topography trend, where k is the number of modes selected during decomposition.

4. The surge correction method based on variational mode decomposition according to claim 1, characterized in that: The seabed reflection wave arrival time curve is subtracted from the seabed topography trend term to obtain the surge correction value on each channel.

5. A surge correction system based on variational mode decomposition, characterized in that: It includes a data acquisition 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 item by variational mode decomposition of the seabed reflection wave arrival time curve; The correction module is used to calculate the surge correction amount on each channel according to the seabed reflection wave arrival time curve and the seabed topography trend item, and perform surge correction using the surge correction amount.

6. A computer storage medium, characterized in that: The computer storage medium stores a computer program, and when the computer program is executed by a processor, the steps of a surge correction method based on variational mode decomposition as claimed in any one of claims 1 to 4 are implemented.

Images