Single-channel seismic data surface multiple suppression method based on wavelet spatial variation compensation

Abstract

The invention relates to the technical field of geophysical exploration signal processing, and relates to a single-channel seismic data surface multiple suppression method based on wavelet space variation compensation. Comprising the steps of constructing a one-dimensional forward modeling model; reconstructing pseudo two-dimensional data and introducing wavelet space variation; constructing an objective function to initialize a primary wave and an operator; carrying out adaptive threshold SVD constraint solution on a surface correlation operator; performing L1 norm constraint inversion to solve primary wave response; judging whether the residual error is converged or not, and if the residual error of the target function is greater than a set value, returning to the optimization algorithm for introducing the L1 norm in the inversion process of the target function; otherwise, jumping out of circulation, and outputting a primary wave result. According to the method, the multiple suppression effect and the stratum interpretation reliability can be remarkably improved, and efficient technical guarantee is provided for safe site selection of clean energy projects such as offshore wind plants; according to the method, the definition and accuracy of the processed seismic section are effectively enhanced, the explanation accuracy is improved, and the reliability of offshore wind plant site assessment is also greatly improved.

Description

Technical Field

[0001] The present invention belongs to the technical field of geophysical exploration signal processing, and in particular relates to a surface multiple wave suppression method for single-channel seismic data based on wavelet spatial variation compensation. Background Art

[0002] With the growing global demand for clean energy, offshore wind power, as a highly efficient form of renewable energy, has developed rapidly in recent years. However, the construction and operation of offshore wind farms face complex geological challenges. In particular, the accurate detection and assessment of seabed geological structures is a key step in ensuring the safety and stability of wind farms. Single-channel seismic (SCS) is well-known for its high efficiency, economy, and high precision, and can be used as the primary means of obtaining seabed geological information. Compared with traditional marine seismic exploration, offshore wind farms are usually located in offshore areas, which results in the development of multiple waves with higher energy in the collected seismic data, causing distortion of the stratum morphology, making it easy for interpreters to give incorrect geological interpretations, and increasing the risk of wind farm site selection.

[0003] Closed-loop surface-related multiple elimination (CL-SRME) is a cutting-edge technology for multiple suppression and a current research hotspot. It directly estimates the primary wave through inversion, eliminating the need for traditional multiple subtraction. This method preserves the integrity of the primary wave signal. Furthermore, its inversion framework is highly stable and has low dependence on the initial model, maintaining robustness in complex noisy environments. Furthermore, it can be applied to data acquired using various seismic acquisition methods, offering broad applicability and flexibility. However, the core algorithm of the traditional CL-SRME method is based on the assumption of spatial invariance of the source wavelet. Since single-channel seismic data only record vertical time series and lack waveform consistency constraints across multiple channels, the wavelet morphology is extremely sensitive to medium heterogeneity. This variability can distort the inversion results. Therefore, it is necessary to account for the spatial variability of the source wavelet. This introduces a more complex inversion problem, making the inversion highly sensitive to noise and unstable, resulting in low inaccuracy and distortion of the final formation location, significantly increasing the risk of offshore wind farm site selection. Summary of the Invention

[0004] The purpose of the present invention is to provide a surface multiple wave suppression method for single-channel seismic data based on wavelet spatial variation compensation, and to solve the problems of low inversion accuracy and poor stability of the traditional CL-SRME method in complex geological environments through the collaborative innovation of wavelet spatial dynamic modeling and adaptive threshold SVD regularization.

[0005] The purpose of the present invention is achieved through the following technical solutions:

[0006] A surface multiple suppression method for single-channel seismic data based on wavelet spatial variation compensation comprises the following steps:

[0007] a. For single-channel seismic data, a one-dimensional primary-multiple forward model is constructed:

[0008]

[0009] Where p(t) represents the seismic acquisition wave field, i.e. the original data, p0(t) represents the primary wave field, represents the multiple wavefield, A represents the surface correlation operator related to the seismic wavelet and reflection coefficient;

[0010] b. Through spatial serial arrangement in the survey line direction, the single-channel data is reconstructed into a pseudo-two-dimensional data profile, providing a data basis for subsequent wavelet spatial variation modeling and multiple wave suppression, that is, P = [p1, p2…, p N ]; Introducing the seismic wavelet A that considers spatial variation into the forward model mult =[A1,A2…,A N ], thus obtaining the primary-multiple forward model considering spatial variation:

[0011] P=P0+P0A mult P (2)

[0012] c. Construct the objective function based on the forward model and obtain the linear relationship between the primary wave and the original data:

[0013]

[0014] Where P and P0 represent the column vector representation of the original data and the primary wave respectively, A represents the surface correlation operator, and I represents the unit operator matrix. Represents the Kronecker product between a matrix and a vector, BlockDiag ω Indicates that a block diagonal matrix is generated on each frequency slice ω, f t and They represent the forward and inverse Fourier transform operators that can realize the mutual transformation between the time domain and the frequency domain, and L is the product of the above-mentioned correlation matrix operators;

[0015] d. Initialize the primary wave response and initialize the surface correlation operator;

[0016] e. In the inversion process of the objective function, the optimization algorithm of the L1 norm is introduced as the constraint condition of the primary wave inversion process to obtain the initial primary wave inversion result;

[0017]

[0018] f. Constructing wavelet matrix model Multi-channel seismic wavelets with spatial correlation are extracted. Their multi-channel joint representation matrix can be regarded as a low-rank structure and is obtained through a low-rank approximation method. An adaptive threshold SVD is introduced to constrain the wavelet model, and the estimated value of the surface correlation operator is solved through the following objective function:

[0019]

[0020] in,‖·‖ * Represents the nuclear norm of the matrix, and uses the nuclear norm of the surface operator as a constraint to fully regularize the solution space;

[0021] g. Substitute the obtained operator estimate into the new primary wave inversion objective function to solve the primary wave response:

[0022]

[0023] Where σ represents the residual between the entire wave field and the estimated full wave field obtained by estimating the first wave;

[0024] h. Calculate the residual of the objective function in step f. If the residual of the objective function is greater than the set value, return to step e. If the residual of the objective function is less than the set value, exit the loop. By alternately iteratively inverting the primary wave and the surface correlation operator, that is, alternately iteratively inverting formula (5) and formula (9), the final primary wave estimation result is directly obtained.

[0025] Furthermore, in step d, the primary wave response is initialized, specifically by applying a cutting function to the original data to obtain an initial primary wave response that does not contain multiple waves.

[0026] Furthermore, in step d, the surface correlation operator is initialized, that is, initially set to a zero matrix.

[0027] Furthermore, in step f, the solution space is fully regularized. The specific solution process is as follows:

[0028] f1, low-rank modeling of the wavelet matrix, whose low-rank waveform structure can be characterized by SVD decomposition:

[0029]

[0030] Among them, U and V represent the left and right singular vector matrices of the wavelet matrix respectively, and the singular value ∑ reflects the energy concentration characteristic of the wavelet matrix;

[0031] f2. Calculate the adaptive threshold according to the wavelet matrix and perform threshold shrinkage:

[0032]

[0033] Among them, λ represents the adaptive threshold, which is obtained through the threshold parameter, and β is the threshold parameter, which is related to the signal-to-noise ratio of the observed data;

[0034] f3, low-rank wavelet estimation update:

[0035]

[0036] Compared with the prior art, the present invention has the following beneficial effects:

[0037] The present invention provides a surface multiple wave suppression method for single-channel seismic data based on wavelet spatial variation compensation. This method can significantly improve the multiple wave suppression effect and the reliability of stratigraphic interpretation, and provides efficient technical guarantees for the safe site selection of clean energy projects such as offshore wind farms. This method effectively enhances the clarity and accuracy of processed seismic profiles, lays a solid foundation for subsequent interpretation, improves the accuracy of interpretation, and greatly improves the reliability of offshore wind farm site assessment. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following briefly introduces the drawings required for use in the embodiments. It should be understood that the following drawings only illustrate certain embodiments of the present invention and therefore should not be regarded as limiting the scope. For ordinary technicians in this field, other relevant drawings can be obtained based on these drawings without paying any creative work.

[0039] Figure 1 This is a flow chart of the steps of the surface multiple wave suppression method of single-channel seismic data based on wavelet spatial variation compensation of the present invention;

[0040] Figure 2 To demonstrate the effect of reconstructing pseudo-two-dimensional data, (a) is single-channel seismic data, (b) is the spatial serialization arrangement, and (c) is the pseudo-two-dimensional seismic data section;

[0041] Figure 3 To demonstrate the resection effect, the initial wave response after resection is performed;

[0042] Figure 4 The following are the changes in the effects of the sub-wave items, where (a) is the change in the effects of the sub-wave items before the constraint items are added, and (b) is the change in the effects of the sub-wave items after the constraint items are added.

[0043] Figure 5 The multiple wave suppression effect is demonstrated, where (a) is the multiple wave suppression effect before the multiple wave suppression, and (b) is the multiple wave suppression effect after the multiple wave suppression. DETAILED DESCRIPTION

[0044] The present invention will be further described below in conjunction with embodiment:

[0045] The present invention will be further described in detail below with reference to the accompanying drawings and examples. It will be understood that the specific embodiments described herein are intended only to illustrate the present invention and are not intended to limit the present invention. It should also be noted that, for ease of description, the accompanying drawings only illustrate portions relevant to the present invention, not all structures.

[0046] It should be noted that similar reference numerals and letters represent similar items in the following drawings. Therefore, once an item is defined in one drawing, it does not need to be further defined or explained in subsequent drawings. At the same time, in the description of the present invention, the terms "first", "second", etc. are used only to distinguish the description and should not be understood as indicating or implying relative importance.

[0047] The present invention improves the traditional CL-SRME method, breaks through the assumption of spatial invariance of the seismic wavelet, constructs a model of dynamic changes of the seismic wavelet with space, and significantly enhances the applicability of the method in complex geological environments. However, it also leads to challenges in the stability of the inversion process. To this end, the present invention innovatively introduces the adaptive threshold SVD (singular value decomposition) technology, integrates it into the inversion process as prior information, and fully regularizes the solution space. This method can effectively capture the intrinsic structural characteristics of seismic data, suppress noise interference through low-rank approximation of the data, and improve the adaptability of the inversion process to the data. Therefore, under complex geological conditions, it not only improves the applicability of the multiple wave suppression method, but also significantly enhances the stability and accuracy of the inversion results.

[0048] The present invention is based on a single-channel seismic data surface multiple wave suppression method based on wavelet spatial variation compensation. Figure 1 As shown, the following steps are included:

[0049] a. For single-channel seismic data, a one-dimensional primary-multiple forward model is constructed:

[0050]

[0051] Where p(t) represents the seismic acquisition wave field, i.e. the original data, p0(t) represents the primary wave field, represents the multiple wavefield, and A represents the surface correlation operator related to the seismic wavelet and reflection coefficient.

[0052] b. Through spatial serial arrangement in the survey line direction, the single-channel data is reconstructed into a pseudo-two-dimensional data profile, providing a data basis for subsequent wavelet spatial variation modeling and multiple wave suppression, that is, P = [p1, p2…, p N ]. Then, the seismic wavelet A considering spatial variation is introduced into the forward model. mult =[A1,A2…,A N], thus obtaining the primary-multiple forward model considering spatial variation:

[0053] P=P0+P0A mult P (2)

[0054] c. Construct the objective function based on the forward model and obtain the linear relationship between the primary wave and the original data:

[0055]

[0056] Where P and P0 represent the column vector representation of the original data and the primary wave respectively, A represents the surface correlation operator, and I represents the unit operator matrix. Represents the Kronecker product between a matrix and a vector, BlockDiag ω Indicates that a block diagonal matrix is generated on each frequency slice ω, f t and They represent the forward and inverse Fourier transform operators that can realize the mutual conversion between time domain and frequency domain.

[0057] d. Initialize the primary wave response (by applying the cutoff function to the original data to obtain the initial primary wave response without multiple waves) and initialize the surface correlation operator (initially set to zero matrix).

[0058] e. In the inversion process of the objective function, the L1 norm optimization algorithm is introduced as the constraint condition of the primary wave inversion process to obtain the initial primary wave inversion result.

[0059]

[0060] f. To effectively extract multi-channel seismic wavelets with spatial correlation, construct a wavelet matrix model Its multi-channel joint representation matrix can be regarded as a low-rank structure and can be obtained by low-rank approximation. The present invention introduces adaptive threshold SVD to constrain the wavelet model and solves the estimated value of the surface correlation operator through the following objective function:

[0061]

[0062] in,‖·‖ * Represents the nuclear norm of the matrix, and uses the nuclear norm of the surface operator as a constraint to fully regularize the solution space. The specific solution process is as follows:

[0063] f1, low-rank modeling of the wavelet matrix, whose low-rank waveform structure can be characterized by SVD decomposition:

[0064]

[0065] The singular value ∑ reflects the energy concentration characteristic of the wavelet matrix.

[0066] f2. Calculate the adaptive threshold according to the wavelet matrix and perform threshold shrinkage:

[0067]

[0068] where β is the threshold parameter, which is related to the signal-to-noise ratio of the observed data.

[0069] f3, low-rank wavelet estimation update:

[0070]

[0071] g. Substitute the obtained operator estimate into the new primary wave inversion objective function to solve the primary wave response:

[0072]

[0073] Where σ is the residual between the entire wavefield and the estimated full wavefield obtained by estimating the first wave.

[0074] h. Calculate the residual of the objective function in step f. If the residual is greater than the set value, return to step e. If the residual is less than the set value, exit the loop. By performing alternating iterative inversion on the primary wave and surface correlation operators, that is, alternating iterative inversion of formula (5) and formula (9), the final primary wave estimation result can be directly obtained.

[0075] Example 1

[0076] The method was applied to the single-channel seismic data measured in the Fangchenggang wind farm in China. The survey line length was 6 km, the channel spacing was 2 m, the sampling interval was 0.005 ms, and the total sampling time was 130 ms.

[0077] This embodiment of the surface multiple suppression method for single-channel seismic data based on wavelet spatial variation compensation is mainly achieved through the following process:

[0078] a. For single-channel seismic data, a one-dimensional primary-multiple wave forward model is constructed.

[0079]

[0080] Where p(t) represents the seismic acquisition wave field, i.e. the original data; p0(t) represents the primary wave field; represents the multiple wavefield, and A represents the surface correlation operator related to the seismic wavelet and reflection coefficient.

[0081] b. Through spatial serial arrangement in the survey line direction, the single-channel data is reconstructed into a pseudo-two-dimensional data profile, providing a data basis for subsequent wavelet spatial variation modeling and multiple wave suppression, that is, P = [p1, p2…, pN ], where p i is a one-dimensional column vector, and the combined P is a two-dimensional matrix. Figure 2 The reconstruction process is shown, in which Figure 2 a is the collected single-channel seismic data, Figure 2 b shows the arrangement of spatial serialization, which intuitively shows the expansion process of spatial serialization on data structure, and obtains pseudo two-dimensional seismic profile ( Figure 2 c), the seismic wavelet has the characteristics of spatial dynamic change. Then, the seismic wavelet A considering spatial change is introduced into the forward model. mult =[A1,A2…,A N ], then the two-dimensional wavelet model A mult Also, the one-dimensional vector A i The permutations and combinations are obtained, thus obtaining the primary-multiple wave forward model considering spatial variations:

[0082] P=P0+P0A mult P (2)

[0083] c. Construct the objective function based on the forward model and obtain the linear relationship between the primary wave and the original data:

[0084]

[0085] Where P and P0 represent the column vector representation of the original data and the primary wave respectively; A represents the surface correlation operator, and I represents the unit operator matrix; Represents the Kronecker product between a matrix and a vector; BlockDiag ω It means that a block diagonal matrix is generated on each frequency slice ω. t and They represent the forward and inverse Fourier transform operators that can realize the mutual conversion between time domain and frequency domain.

[0086] d. Initialize the wave response by applying the cut function to the original data ( Figure 2 c) Obtain the initial primary wave response without multiple waves and initialize the surface correlation operator (initial setting is zero matrix). Figure 3 As shown in FIG, in this experiment, the data outside the sampling points 920-1010 are cut off, and only the primary wave energy received before the multiple waves appear is retained within the sampling points, so as to obtain the initial data containing only the primary wave response.

[0087] e. In the inversion process of the objective function, the L1 norm optimization algorithm is introduced as the constraint condition of the primary wave inversion process to obtain the initial primary wave inversion result.

[0088] In this embodiment, the L1-norm optimization algorithm is reflected in the constraint term of formula (4), namely, ||P0||1≤τ. Its purpose is to constrain the objective function by calculating the L-norm of the initial primary wave, thereby retaining characteristic information, removing redundant information, and making the inversion process more stable. Therefore, the L1-norm constraint optimization process can be reflected in formula (4).

[0089]

[0090] In this experiment, τ is obtained based on the L1 norm of the initial primary wave, which is 30.358.

[0091] f. To effectively extract multi-channel seismic wavelets with spatial correlation, construct a wavelet matrix model Its multi-channel joint representation matrix can be regarded as a low-rank structure and can be obtained by a low-rank approximation method. This embodiment introduces an adaptive threshold SVD to constrain the wavelet model and obtains the estimated value of the surface correlation operator through the following objective function:

[0092]

[0093] in,‖·‖ * Represents the nuclear norm of the matrix, and uses the nuclear norm of the surface operator as a constraint to fully regularize the solution space. The specific solution process is as follows:

[0094] f1, low-rank modeling of the wavelet matrix, whose low-rank waveform structure can be characterized by SVD decomposition:

[0095]

[0096] The singular value ∑ reflects the energy concentration characteristic of the wavelet matrix.

[0097] f2. Calculate the adaptive threshold according to the wavelet matrix and perform threshold shrinkage:

[0098]

[0099] Where β is a threshold parameter, which is related to the signal-to-noise ratio of the observed data. In this experiment, β = 0.15.

[0100] f3, low-rank wavelet estimation update:

[0101]

[0102] Figure 4 The change of the wavelet morphology after adding regularization constraints is shown. It can be seen that the existence of bad tracks before the constraints may affect the results of the wavelet terms, resulting in poor lateral continuity of the wavelet terms, such as Figure 4The constrained wavelet morphology significantly enhances the spatial lateral continuity and eliminates the impact of individual bad tracks in the data on the suppression of multiple waves, as shown in Figure 2. Figure 4 As shown in b.

[0103] g. Substitute the obtained operator estimate into the new primary wave inversion objective function to solve the primary wave response:

[0104]

[0105] Wherein, the value of σ is 1%-10% of the L2 norm of the original data, and its physical meaning is the residual between the entire wave field and the estimated full wave field obtained by the estimated primary wave calculation.

[0106] h. Calculate the residual of the objective function in step f. If the residual of the objective function is greater than the set value, return to step e; if the residual of the objective function is less than the set value, the loop can be jumped out. By alternately iteratively inverting the primary wave and the surface correlation operator, formula (5) and formula (9) are alternately iteratively inverted to directly obtain the final primary wave estimation result. The processing result is as follows Figure 5 As shown, it can be seen that in Figure 5 In the black frame in a, there are also false strata composed of multiple waves. Figure 5 Figure b shows the results of this method, where the false strata originally composed of multiple waves are suppressed. Furthermore, because noise is also suppressed to a certain extent, the original strata are easier to identify, preventing subsequent interpretation personnel from misidentifying the strata and potentially causing construction risks.

[0107] Note that the above are only preferred embodiments of the present invention and the technical principles employed. Those skilled in the art will understand that the present invention is not limited to the specific embodiments described herein, and that various obvious changes, readjustments, and substitutions can be made by those skilled in the art without departing from the scope of protection of the present invention. Therefore, although the present invention has been described in detail through the above embodiments, the present invention is not limited to the above embodiments and may include many other equivalent embodiments without departing from the concept of the present invention. The scope of the present invention is determined by the scope of the appended claims.

Claims

1. A surface multiple suppression method for single-channel seismic data based on wavelet spatial variation compensation, characterized in that: The following steps are involved: a. For single-channel seismic data, a one-dimensional primary-multiple forward model is constructed: Where p(t) represents the seismic acquisition wave field, i.e. the original data, p0(t) represents the primary wave field, represents the multiple wavefield, A represents the surface correlation operator related to the seismic wavelet and reflection coefficient; b. Through spatial serial arrangement in the survey line direction, the single-channel data is reconstructed into a pseudo-two-dimensional data profile, providing a data basis for subsequent wavelet spatial variation modeling and multiple wave suppression, that is, P = [p1, p2…, p N ]; Introducing the seismic wavelet A that considers spatial variation into the forward model mult =[A1,A2…,A N ], thus obtaining the primary-multiple forward model considering spatial variation: P=P0+P0A mult P (2) c. Construct the objective function based on the forward model and obtain the linear relationship between the primary wave and the original data: Where P and P0 represent the column vector representation of the original data and the primary wave respectively, A represents the surface correlation operator, and I represents the unit operator matrix. Represents the Kronecker product between a matrix and a vector, BlockDiag ω Indicates that a block diagonal matrix is generated on each frequency slice ω, f t and They represent the forward and inverse Fourier transform operators that can realize the mutual transformation between the time domain and the frequency domain respectively; d. Initialize the primary wave response and initialize the surface correlation operator; e. In the inversion process of the objective function, the L1 norm optimization algorithm is introduced as the constraint condition of the primary wave inversion process to obtain the initial primary wave inversion result. f. Constructing wavelet matrix model Multi-channel seismic wavelets with spatial correlation are extracted. Their multi-channel joint representation matrix can be regarded as a low-rank structure and is obtained through a low-rank approximation method. An adaptive threshold SVD is introduced to constrain the wavelet model, and the estimated value of the surface correlation operator is solved through the following objective function: in,‖·‖ * Represents the nuclear norm of the matrix, and uses the nuclear norm of the surface operator as a constraint to fully regularize the solution space; g. Substitute the obtained operator estimate into the new primary wave inversion objective function to solve the primary wave response: Where σ represents the residual between the entire wave field and the estimated full wave field obtained by estimating the first wave; h. Calculate the residual of the objective function in step f. If the residual of the objective function is greater than the set value, return to step e. If the residual of the objective function is less than the set value, exit the loop. By alternately iteratively inverting the primary wave and the surface correlation operator, that is, alternately iteratively inverting formula (5) and formula (9), the final primary wave estimation result is directly obtained.

2. The surface multiple suppression method for single-channel seismic data based on wavelet spatial variation compensation according to claim 1, characterized in that: Step d: Initializing the primary wave response, specifically, applying a cutting function to the original data to obtain an initial primary wave response that does not contain multiple waves.

3. The surface multiple suppression method for single-channel seismic data based on wavelet spatial variation compensation according to claim 1, characterized in that: Step d: Initialize the surface correlation operator, that is, initially set it to a zero matrix.

4. The surface multiple suppression method for single-channel seismic data based on wavelet spatial variation compensation according to claim 1, characterized in that: Step f: Fully regularize the solution space. The specific solution process is as follows: f1, low-rank modeling of the wavelet matrix, whose low-rank waveform structure can be characterized by SVD decomposition: The singular value ∑ reflects the energy concentration characteristic of the wavelet matrix; f2. Calculate the adaptive threshold according to the wavelet matrix and perform threshold shrinkage: Among them, β is the threshold parameter, which is related to the signal-to-noise ratio of the observed data; f3, low-rank wavelet estimation update:

Images