A multi-channel bipolar sparse spectral inversion method based on hessian matrix constraint

Abstract

The present application belongs to the field of oil and gas and geophysical exploration, and particularly relates to a multi-channel bipolar sparse spectral inversion method based on Hessian matrix constraint. The method comprises the following steps: estimating a seismic wavelet wavelength to determine a time window length, and dividing multi-channel seismic data; extracting a seismic wavelet in the time window; constructing a multi-channel bipolar sparse spectral inversion objective function based on Hessian matrix constraint in the time window, and calculating the reflection coefficient of the multi-channel seismic data in the time window; and combining the reflection coefficient results of each time window center in time sequence to obtain the spectral inversion result of the multi-channel seismic data. The method can reflect the real contact relationship of the stratum, depict the real pinch-out position of the stratum, and provide high-resolution and more abundant geological information. While improving the resolution of the seismic data, the method can ensure the lateral continuity of the stratum in the multi-channel data, and can provide more accurate and abundant geological and reservoir information.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of oil and gas and geophysical exploration, and particularly relates to a multi-channel bipolar sparse spectral inversion method based on Hessian matrix constraint. BACKGROUND

[0002] Due to the influence of tuning effect, conventional seismic data and traditional processing methods cannot identify thin layers and small-scale geological structures with a thickness less than the tuning thickness, and it is difficult to provide high-precision interpretation results. The spectral inversion method does not rely on prior information such as logging data, and only needs to use the spectral information of the seismic data itself to invert the sparse reflection coefficient of the stratum. Since it can distinguish seismic thin layers below the tuning thickness, the method has been widely studied and discussed.

[0003] The traditional spectral inversion method uses a single seismic trace to globally calculate the spectrum, and adds a 1-norm constraint term to the objective function for solving. Calculating the spectrum globally for a single seismic trace may cause partial frequency domain signals to overlap, making it impossible to highlight local effective information, and ignoring the non-stationarity in the seismic wave propagation process, which is likely to cause errors.

[0004] Adding a 1-norm constraint term to the objective function of a single seismic trace can sparsely constrain the reflection coefficient sequence, but since multiple seismic traces are calculated separately, the number of iterations for algorithm convergence is different, and the reflection coefficients of each trace have no spatial continuity, resulting in phenomena such as "misplacement" and "bifurcation" on the reflection coefficient profile, which cannot reflect the true distribution of the stratum reflection coefficient.

[0005] Chinese Patent Application CN111856559A discloses a multi-channel seismic spectral inversion method based on sparse Bayesian learning theory, which includes: obtaining preprocessed post-stack seismic data; extracting multi-channel seismic spectral inversion parameters from the post-stack seismic data; extracting a seismic wavelet from the post-stack seismic data; extracting spectral information of the seismic wavelet based on Fourier transform; constructing a forward operator matrix of multi-channel seismic spectral inversion according to the multi-channel seismic spectral inversion parameters; extracting multi-channel average seismic trace spectral information from the post-stack seismic data based on Fourier transform; constructing a multi-channel seismic spectral inversion objective function based on the spectral information, and solving the objective function by sparse Bayesian learning theory based on the maximum expectation algorithm to obtain a multi-channel seismic spectral inversion result.

[0006] Chinese patent application CN111796325A discloses a frequency division iteration constrained stochastic inversion method, including: selecting a drilled well needed for inversion; fine well-seismic calibration; expanding the existing seismic data frequency band; performing frequency division processing on the wide frequency seismic data after expansion; extracting high frequency wavelet; obtaining a three-dimensional apparent sandstone-to-ground ratio distribution body; a seismic data driven inversion method; improving the accuracy of the seismic data driven stochastic inversion result, adding the obtained three-dimensional apparent sandstone-to-ground ratio distribution body as prior information into the high frequency seismic data driven stochastic inversion process, thereby obtaining a high-precision frequency division iteration constrained stochastic inversion result; the method effectively improves the reservoir prediction accuracy, more accurately depicts the lateral thickness variation of the reservoir, restores the real superimposed connected relationship of the reservoir, and provides an important reference basis for pre-drilling well site deployment in the exploration target evaluation stage and well site optimization in the oilfield comprehensive adjustment stage.

[0007] However, based on the fact that each trace reflection coefficient has no spatial coherence, the phenomenon of "misplacement" and "bifurcation" occurs on the reflection coefficient profile, which cannot reflect the real distribution of the formation reflection coefficient. SUMMARY

[0008] The main purpose of the present application is to provide a multi-channel bipolar sparse spectral inversion method based on Hessian matrix constraint, which can improve the resolution of seismic data while ensuring the lateral continuity of the formation in multi-channel data, and can provide more accurate and rich geological and reservoir information, overcoming the shortcomings of the prior art.

[0009] To achieve the above purpose, the technical scheme adopted by the present application is as follows:

[0010] The present application provides a multi-channel bipolar sparse spectral inversion method based on Hessian matrix constraint, comprising the following steps:

[0011] Estimating the wavelength of the seismic wavelet to determine the time window length, and dividing the multi-channel seismic data; extracting the seismic wavelet in the time window; constructing a multi-channel bipolar sparse spectral inversion objective function based on Hessian matrix constraint in the time window, and calculating the reflection coefficient of the multi-channel seismic data in the time window; combining the reflection coefficient results of each time window center in time sequence to obtain the spectral inversion result of the multi-channel seismic data.

[0012] Further, the specific method for estimating the wavelength of the seismic wavelet to determine the time window length is as follows:

[0013] Estimating the main frequency f of the wavelet of the multi-channel seismic data s(x, t) m :

[0014]

[0015] where S(x,f) represents the amplitude spectrum of the multi-channel signal s(x,t), and M is the number of seismic channels;

[0016] The main frequency is f w The length of the zero-phase Ricker wavelet is used as the window length of the time window, so that the sampling point at the center of the time window corresponds to the main peak position of the zero-phase Ricker wavelet.

[0017] Further, the method for dividing the multi-channel seismic data is: taking a time sampling interval as a window moving step, and using the time window to divide the expanded seismic data of each channel; the center point of each time window corresponds to the sampling point of the original data one by one.

[0018] Further, the statistical wavelet estimation method is used in the time window to extract the seismic wavelet with time-varying characteristics:

[0019] w(t)=[1-2(πf m t) 2 ]exp[-(πf m t) 2 ]

[0020] where f m is the main frequency of the Ricker wavelet estimated in the time window.

[0021] Further, a multi-channel bipolar sparse spectrum inversion objective function based on Hessian matrix constraint is constructed:

[0022] The convolution operation of the multi-channel seismic convolution model in the time domain can be transformed into the multiplication operation in the frequency domain:

[0023]

[0024] The objective function of the spectrum inversion can be constructed in a single time window:

[0025] obj=W(f)×R(x,f)-S(x,f)

[0026]

[0027] where,

[0028] A is a matrix constructed by the wavelet spectrum through odd-even decomposition:

[0029]

[0030] A 11 ={α o Re[W(f i )]sin(m i,k )sin(n i )-α o Im[W(f i )]sin(mi,k )cos(n i )}

[0031] A 12 ={α e Re[W(f i )]cos(m i,k )cos(n i )+α e Im[W(f i )]cos(m i,k )sin(n i )}

[0032] A 21 ={α o Re[W(f i )]sin(m i,k )cos(n i )+α o Im[W(f i )]sin(m i,k )sin(n i )}

[0033] A 22 ={-α e Re[W(f i )]cos(m i,k )sin(n i )+α e Im[W(f i )]cos(m i,k )cos(n i )}

[0034] m ij =πf i T j

[0035] n i =πf i (N-1)Δt

[0036] b is the matrix constructed by odd-even decomposition of the spectrum of multi-channel seismic records:

[0037]

[0038] r o 、r e are matrices composed of odd and even components of multi-channel reflection coefficients respectively. By merging odd and even components, the multi-channel reflection coefficient matrix r(x,t) can be obtained:

[0039]

[0040] α o, α e are odd and even weighting coefficients respectively, N is the number of sampling points in the time window, M is the number of channels, H(r) is the Hessian matrix of the reflection coefficient r, ||g|| 2，1 is the L2,1 norm; ||g|| F is the Frobenius or Euclidean norm.

[0041] Furthermore, an optimal solution algorithm is used to find the minimum value of the spectral inversion objective function within the time window. The objective function of the multi-channel bipolar sparse spectral inversion based on Hessian matrix constraints is a typical linear regression equation. Accurate spectral inversion results can be obtained through iterative fitting using conventional least squares methods. Global optimization algorithms such as genetic algorithms, simulated annealing, and particle swarm optimization can also yield accurate spectral inversion results.

[0042] Furthermore, all the segmented time windows must be calculated.

[0043] Furthermore, the reflection coefficient values ​​of the center points of each time window are sorted along the time sequence, and finally the ideal sparse spectrum inversion result is obtained:

[0044] r(x,t)=[r1(x,t)·δ(x,t) L r P (x,t)·δ(x,t)] M×P

[0045] where r P (x, t) is the multichannel reflection coefficient calculated by spectral inversion within the Pth time window, δ(x, t) is the Dirac function matrix, M is the number of channels, and P is the number of time windows. Recording the values ​​of the resulting multichannel reflection coefficient sequence at the central sampling point of the time window—that is, applying a Dirac function window to the multichannel reflection coefficient result from spectral inversion within the time window—can eliminate the effects of spectral leakage caused by time window segmentation.

[0046] Compared with the prior art, the present invention has the following advantages:

[0047] The proposed method adds time windows to multiple seismic channels, extracts local signals, and calculates the spectrum simultaneously. The signals within the time window are approximately stationary, highlighting their local spectral characteristics. The extracted wavelets exhibit both spatially and time-varying characteristics. Simultaneous calculation of multiple channels ensures parameter consistency across all channels.

[0048] The method of the present invention also adds a Hessian matrix of the reflection coefficient to the spectral inversion objective function to perform structural constraints. The Hessian matrix elements are composed of the second-order partial derivatives of the reflection coefficient with respect to space, which can be used to detect three-dimensional spatial interface singularities. The Frobenius or Euclidean norm of the Hessian matrix can determine the hypersurface corresponding to the layer boundary and fault. The Hessian matrix constraint forces the inverted reflection coefficient to be distributed along a smooth surface in 3D space, that is, the reflection coefficient has good lateral continuity between each channel, which is more consistent with the distribution characteristics of the real formation reflection coefficient.

[0049] The method of the present invention can separate the complex wave events formed by coherence, identifying thin layers and small-scale geological bodies whose thickness is less than the tuning thickness that are indistinguishable on the original seismic profile. This method reflects the true contact relationship of the strata and depicts the true pinch-out location of the strata, providing high-resolution and richer geological information. While improving the resolution of seismic data, this method can also ensure the lateral continuity of strata in multi-channel data, providing more accurate and richer geological and reservoir information; therefore, the method of the present invention can play an important role in reservoir prediction and oil and gas detection. BRIEF DESCRIPTION OF THE DRAWINGS

[0050] The accompanying drawings, which constitute a part of the present invention, are used to provide a further understanding of the present invention. The exemplary embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute improper limitations on the present invention.

[0051] Figure 1 This is a flow chart of a multi-channel bipolar sparse spectrum inversion method based on Hessian matrix constraints according to a specific embodiment of the present invention;

[0052] Figure 2 This is a seismic forward section of the sand-mud interbedded wedge model according to a specific embodiment of the present invention;

[0053] Figure 3 A reflection coefficient profile obtained by using a multi-channel bipolar sparse spectrum inversion method based on Hessian matrix constraints for the seismic forward modeling profile of the sand-mud interbedded wedge model described in a specific embodiment of the present invention;

[0054] Figure 4 A reflection coefficient profile obtained by using a conventional spectral inversion method for a seismic forward model of a sand-mud interbedded wedge model according to a specific embodiment of the present invention;

[0055] Figure 5 A seismic forward modeling section of the Marmousi model according to a specific embodiment of the present invention;

[0056] Figure 6The reflection coefficient profile used in the Marmousi model seismic forward profile of one embodiment of the present application is obtained by using a multi-channel bipolar sparse spectral inversion method based on Hessian matrix constraint.

[0057] Figure 7 The reflection coefficient profile used in the Marmousi model seismic forward profile of one embodiment of the present application is obtained by using a multi-channel bipolar sparse spectral inversion method based on Hessian matrix constraint. DETAILED DESCRIPTION

[0058] It should be noted that the following detailed description is exemplary in nature and is intended to provide further description of the application. Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs.

[0059] It is to be understood that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting of example embodiments in accordance with the present application. As used herein, the singular forms "a", "an" and "the" are intended to include the plural forms as well, unless the context clearly indicates otherwise. It will be further understood that the terms "comprises" and / or "comprising," when used in this specification, specify the presence of stated features, steps, operations, elements, and / or components, but do not preclude the presence or addition of one or more other features, steps, operations, elements, components, and / or groups thereof.

[0060] In order to enable a person skilled in the art to more clearly understand the technical solutions of the present application, the technical solutions of the present application will be described in detail below with specific embodiments.

[0061] Embodiment 1

[0062] As shown in the following table, the multi-channel bipolar sparse spectral inversion method based on Hessian matrix constraint includes the following steps: Figure 1

[0063] S1, input multi-channel seismic data.

[0064] S2, estimate the seismic wavelet to determine the window length of the time window, and divide the data by using the time window:

[0065] First, estimate the wavelet main frequency f m :

[0066]

[0067] Where S(x,f) represents the amplitude spectrum of the multi-channel signal s(x,t), and M is the number of seismic channels.

[0068] Take the main frequency f w The zero-phase Ricker wavelet length is taken as the window length of the time window, so that the center sampling point of the time window corresponds to the main peak position of the zero-phase Ricker wavelet.

[0069] ​After that, the window step size is set to a time sampling interval of the multi-channel seismic data, and each channel of seismic data is divided using the time window. At this time, the center point of each time window corresponds to each sampling point of the original data.

[0070] S3. Use statistical wavelet estimation method within the time window to extract seismic wavelets with time-varying characteristics:

[0071] w(t)=[1-2(πf m t) 2 ]exp[-(πf m t) 2 ],

[0072] Among them, f m is the estimated dominant frequency of the Ricker wavelet within the time window.

[0073] S4. Construct and solve the multi-channel bipolar sparse spectrum inversion objective function based on Hessian matrix constraints:

[0074] The convolution operation of the multi-channel seismic convolution model in the time domain can be transformed into a product operation in the frequency domain.

[0075]

[0076] The objective function of spectral inversion can be constructed within a single time window:

[0077] obj=W(f)×R(x,f)-S(x,f)

[0078]

[0079] in,

[0080] A is the matrix constructed by odd-even decomposition of the wavelet spectrum:

[0081]

[0082] A 11 ={α o Re[W(f i )]sin(m i,k )sin(n i )-α o Im[W(f i )]sin(m i,k )cos(n i )}

[0083] A 12 ={α e Re[W(f i )]cos(m i,k )cos(n i )+αe Im[W(f i )]cos(m i,k )sin(n i )}

[0084] A 21 ={α o Re[W(f i )]sin(m i,k )cos(n i )+α o Im[W(f i )]sin(m i,k )sin(n i )}

[0085] A 22 ={-α e Re[W(f i )]cos(m i,k )sin(n i )+α e Im[W(f i )]cos(m i,k )cos(n i )}

[0086] m ij =πf i T j

[0087] n i =πf i (N-1)Δt

[0088] b is the matrix constructed by odd-even decomposition of the spectrum of multi-channel seismic records:

[0089]

[0090] r o 、r e are matrices composed of odd and even components of multi-channel reflection coefficients respectively. By merging odd and even components, the multi-channel reflection coefficient matrix r(x,t) can be obtained:

[0091]

[0092] α o , α e are odd and even weighting coefficients respectively, N is the number of sampling points in the time window, M is the number of channels, H(r) is the Hessian matrix of the reflection coefficient r, ||g|| 2，1 is the L2,1 norm; ||g|| F is the Frobenius or Euclidean norm.

[0093] The objective function of the multi-channel bipolar sparse spectral inversion based on Hessian matrix constraints is a typical linear regression equation. Accurate spectral inversion results can be obtained by iteratively fitting and solving the minimum value of the spectral inversion objective function using conventional least-squares methods. Global optimization algorithms such as genetic algorithms, simulated annealing, and particle swarm optimization can also yield accurate spectral inversion results.

[0094] S5. Repeat steps S3 and S4 until all time windows divided in step S2 are calculated.

[0095] S6. The reflection coefficient values ​​of the center points of each time window are calculated in time sequence to obtain the ideal sparse spectrum inversion result: r(x,t)=[r1(x,t)·δ(x,t) L r P (x,t)·δ(x,t)] M×P

[0096] where r P (x, t) is the multichannel reflection coefficient calculated by spectral inversion within the Pth time window, δ(x, t) is the Dirac function matrix, M is the number of channels, and P is the number of time windows. Recording the values ​​of the resulting multichannel reflection coefficient sequence at the central sampling point of the time window—that is, applying a Dirac function window to the multichannel reflection coefficient result from spectral inversion within the time window—can eliminate the effects of spectral leakage caused by time window segmentation.

[0097] Example 2

[0098] The seismic forward section of the sand-mud interbedded wedge model is inverted using the multi-channel bipolar sparse spectrum inversion method based on Hessian matrix constraints described in Example 1. The seismic forward section of the sand-mud interbedded wedge model is as follows: Figure 2 As shown, the reflection coefficient profile is calculated as Figure 3 shown.

[0099] Example 3

[0100] The Marmousi model seismic forward section is inverted using the multi-channel bipolar sparse spectrum inversion method based on Hessian matrix constraints described in Example 1. The Marmousi model seismic forward section is as follows: Figure 5 As shown, the reflection coefficient profile is calculated as Figure 6 shown.

[0101] Will Figure 3 、 Figure 6 The inversion results shown are consistent with the conventional spectrum inversion results ( Figure 4 、 Figure 7Compared with the conventional spectral inversion, the thin layer resolution capability is relatively poor, the thin layer with a thickness less than 1 / 8 wavelength cannot be distinguished, and the reflection coefficient position is shifted and lacks lateral continuity; the multi-channel bipolar sparse spectral inversion based on the Hessian matrix constraint can distinguish the thin layer with a thickness less than 1 / 8 wavelength, and the reflection coefficient has good lateral continuity between channels; the spectral inversion result of the application has higher thin layer resolution capability, is more in line with the distribution characteristics of the real formation reflection coefficient, and has high application value.

[0102] The above embodiment is a preferred embodiment of the present application, but the embodiment of the present application is not limited by the above embodiment, and any change, modification, substitution, combination, simplification made without departing from the spirit and principle of the present application should be an equivalent replacement mode, and all are included in the protection scope of the present application.

Claims

1. A multi-channel bipolar sparse spectrum inversion method based on Hessian matrix constraints, characterized in that: The following steps are involved: Estimate the wavelength of the seismic wavelet to determine the time window length and segment the multi-channel seismic data; extract the seismic wavelet within the time window; construct a multi-channel bipolar sparse spectral inversion objective function based on the Hessian matrix constraint within the time window and calculate the reflection coefficient of the multi-channel seismic data within the time window; combine the reflection coefficient results of the center of each time window along the time sequence to obtain the spectral inversion result of the multi-channel seismic data; The objective function of spectral inversion is constructed within a single time window: obj=W(f)×R(x,f)-S(x,f) in, is the dipole spectrum inversion loss function, which describes the distribution characteristics of the seismic data noise spectrum using the reflection coefficient of odd-even decomposition based on Bayesian theory; λ||[r(x,t)]|| 2，1 It is a multi-channel sparse constraint term for dipole spectrum inversion, which effectively describes the three-dimensional group sparse characteristics of multi-channel reflection coefficients, that is, the distribution characteristics of horizontal continuity and vertical sparseness; μ||H[r(x,t)]|| F It is the Hessian matrix constraint term. The hypersurface of the real formation reflection interface can be determined through the Frobenius norm of the Hessian matrix, so that the reflection coefficient result based on the Hessian matrix constraint inversion is more consistent with the distribution characteristics of the real formation reflection coefficient. A is the matrix constructed by odd-even decomposition of the wavelet spectrum: A 11 ={a o Re[W(f i )]sin(m i,j )sin(n i )-a o Im[W(f i )]sin(m i,j )cos(n i )} A 12 ={a e Re[W(f i )]cos(m i,j )cos(n i )+a e Im[W(f i )]cos(m i,j )sin(n i )} A 21 ={a o Re[W(f i )]sin(m i,j )cos(n i )+a o Im[W(f i )]sin(m i,j )sin(n i )} A 22 ={-a e Re[W(f i )]cos(m i,j )sin(n i )+a e Im[W(f i )]cos(m i,j )cos(n i )} m ij =πf i T j n i =πf i (N-1)Δt b is the matrix constructed by odd-even decomposition of the spectrum of multi-channel seismic records: α o , α e are odd and even weighted coefficients respectively, N is the number of sampling points in the time window, M is the number of channels, H(r) is the Hessian matrix of the reflection coefficient r, || || 2.1 is the L2,1 norm; || || F is the Frobenius or Euclidean norm.

2. The method according to claim 1, characterized in that Specific method for estimating the wavelength of seismic wavelets and determining the time window length: Estimate the dominant frequency f of the wavelet of multi-channel seismic data s(x,t) m : Where S(x,f) represents the amplitude spectrum of the multi-channel signal s(x,t), and M is the number of seismic channels; Take the main frequency as f w The length of the zero-phase Ricker wavelet is used as the window length of the time window, so that the central sampling point of the time window corresponds to the main peak position of the zero-phase Ricker wavelet.

3. The method according to claim 1, characterized in that The method for segmenting multi-channel seismic data is as follows: using a time sampling interval as the window step, and using the time window to segment the expanded seismic data; the center point of each time window corresponds to the sampling point of the original data one by one.

4. The method according to claim 1, characterized in that The statistical wavelet estimation method is used within the time window to extract the seismic wavelet with time-varying characteristics: w(t)=[1-2(πf m t) 2 ]exp[-(πf m t) 2 ] Among them, f m is the estimated dominant frequency of the Ricker wavelet within the time window.

5. The method according to claim 1, characterized in that: The optimal solution algorithm is used to find the minimum value of the spectral inversion objective function within the time window.

6. The method according to claim 1, characterized in that All the split time windows must be calculated.

7. The method according to claim 1, characterized in that: The reflection coefficient values ​​of the center points of each time window are calculated in time sequence to obtain the ideal sparse spectrum inversion result: r(x,t)=[r1(x,t)·δ(x,t) … r P (x,t)·δ(x,t)] M×P where r P (x, t) is the multi-channel reflection coefficient calculated by spectrum inversion in the Pth time window, δ(x, t) is the Dirac function matrix, M is the number of channels, and P is the number of time windows; Recording the values ​​of the obtained multi-channel reflection coefficient sequence at the sampling point in the center of the time window, that is, adding a Dirac function window to the multi-channel reflection coefficient result of the spectrum inversion within the time window, can eliminate the spectrum leakage effect caused by time window division.

Images