Blind inversion method based on combined constraint of construction and non-convex total variation

Abstract

The invention provides a blind inversion method based on construction and non-convex total variation combined constraint, and belongs to the field of seismic inversion. The method comprises the following steps: inputting seismic data, and determining a to-be-inverted post-stack seismic section; calculating an initial seismic wavelet and an initial wave impedance model; calculating a seismic dip angle through the post-stack seismic data, and constructing a rotation operator; an objective function of wave impedance inversion is established, a split-bregman algorithm is adopted for solving, and wave impedance data are output; an objective function of wavelet inversion is established, wave impedance data are substituted, a split-bregman algorithm is adopted for solving, and seismic wavelets are output; and alternately updating the target functions of the wave impedance and the wavelets, and when an iteration criterion is satisfied, outputting the inverted wave impedance and the seismic wavelets. The invention provides a frame for simultaneously inverting the wave impedance and the wavelets, and the introduced structural constraint and non-convex total variation regularization can effectively improve the accuracy of the wave impedance and wavelet inversion.

Description

Technical Field

[0001] The invention relates to the field of seismic exploration of coal fields, oil and gas and other mineral resources, and in particular to a blind inversion method based on joint constraints of structure and non-convex total variation. Background Art

[0002] Seismic exploration is an important method for coalfield and oil and gas resource exploration. Wave impedance inversion is the main means to obtain underground wave impedance information through seismic data. It converts seismic data into a model space with higher resolution, so that an accurate wave impedance model can be constructed. Conventional wave impedance inversion methods use total variation (TV) regularization to invert underground wave impedance. TV regularization uses L1 norm constraints on the horizontal and vertical derivatives of model parameters to construct the objective function. Zhang (2014) et al. used L1 norm mismatch function, TV regularization and prior information constraints to establish seismic inversion objective function through Lagrange multiplier method, solving the inversion problem of outliers and stratum boundary discontinuities in seismic data. Gholami (2015) used TV regularization to recover the block wave impedance model from reflectivity data and improved the spatiotemporal correlation of model parameters. Wang et al. (2018) proposed that compared with isotropic total variation (ITV), anisotropic total variation (ATV) better preserves the discontinuous interface between strata. Zhang et al. (2022) developed an impedance inversion method based on geological structure-guided TV (GSGTV) constraints with improved resolution and spatial continuity.

[0003] However, seismic wavelets are unknown. Generally, theoretical Ricker wavelets or statistical wavelets can be used for inversion, but this will bring certain errors. Therefore, Gholami (2016) developed a blind inversion method to simultaneously recover wavelet and wave impedance information from seismic data. However, TV regularization cannot fully exploit sparsity, has inherent bias, and cannot accurately recover stratum boundary information. In addition, TV regularization does not take into account geological structure and is not suitable for seismic inversion in complex areas. Summary of the invention

[0004] In view of the shortcomings of existing blind inversion methods, the present invention proposes a blind inversion method based on joint constraints of construction and non-convex total variation, which mainly solves the problems of inherent bias in TV regularization, inability to accurately restore stratum boundary information and unsuitability for complex areas.

[0005] In order to solve the above technical problems, the technical solution adopted by the present invention is: a blind inversion method based on construction and non-convex total variation joint constraints, comprising:

[0006] Step 1: Input seismic data and calculate the initial seismic wavelet and initial wave impedance model;

[0007] Step 2: Calculate the seismic dip angle through post-stack seismic data and construct a rotation operator;

[0008] Step 3: Establish the objective function of wave impedance inversion, bring in the initial seismic wavelet, and use the split-bregman algorithm to solve and obtain the wave impedance model;

[0009] Step 4: Establish the objective function of wavelet inversion, bring in the wave impedance model, and use the split-bregman algorithm to solve and obtain the seismic wavelet;

[0010] Step 5: Determine whether the iteration criteria are met. If not, return to step 3. If so, output the final wave impedance and wavelet inversion results.

[0011] Step 1 includes the following:

[0012] Input seismic data, and perform stacking processing on the seismic data to obtain post-stack seismic data;

[0013] Perform spectrum analysis on the seismic data, estimate the main frequency of the seismic data, and calculate the zero-phase Ricker wavelet of the corresponding main frequency as the initial seismic wavelet:

[0014] The initial seismic wavelet is convolved with the reflection coefficient sequence calculated from the well logging curve to obtain a synthetic seismic record, which is then compared with the actual seismic record until the main stratigraphic positions match, and then fine-tuned to maximize the correlation coefficient, thereby obtaining the well logging curve in the time domain;

[0015] Starting from the logging curve, the stratigraphic layer is tracked along the seismic phase axis, and a two-dimensional profile is obtained by interpolating the one-dimensional logging data along the layer. The interpolated logarithmic wave impedance is low-pass filtered with a cutoff frequency of 5 Hz to obtain the initial wave impedance model, and its logarithm is calculated.

[0016] Step 2 includes the following:

[0017] The gradient structure tensor is calculated based on the seismic data, the local dip is calculated from the gradient structure tensor, and the rotation operator is constructed from the local dip:

[0018]

[0019] D perp is a first-order difference rotation operator perpendicular to the local structure direction, D parl is a first-order difference rotation operator parallel to the local structure. cos and Q sin are diagonal matrices consisting of the cosine and sine of the local inclination angle, respectively.

[0020] Step 3 includes the following:

[0021] Establish the objective function of wave impedance inversion:

[0022]

[0023] subjecttoR perp =D perp Z,R parl =D parl Z

[0024] Where G = WD, W is the wavelet matrix, D is the 0.5 times first-order difference matrix, Z is the logarithmic impedance to be inverted, Φ α () is the non-convex smooth clipped absolute deviation (SCAD) penalty, R perp and R parl are auxiliary variables, μ and σ are regularization parameters;

[0025] The split-bregman algorithm is used to solve the objective function and obtain the wave impedance model.

[0026] Step 4 includes the following:

[0027] Establish the objective function of wavelet inversion:

[0028]

[0029] in is the seismic sub-wave vector, R is the diagonal matrix composed of the Toeplitz matrix of each reflection coefficient sequence, is the initial seismic wavelet, D1 and D2 are the first-order difference matrix and the second-order difference matrix respectively. μ1, μ2, μ3, and μ4 are regularization parameters.

[0030] The split-bregman algorithm is used to solve the objective function and obtain the seismic wavelet.

[0031] Step 5 includes the following:

[0032] Define the iteration stopping criterion:

[0033]

[0034] The superscript k+1 indicates the result of the k+1th cycle in the wavelet inversion cycle;

[0035] Given the iteration stop parameter tol3, determine whether the external iteration stop criterion is met: if tol3 is less than TOL3, then let Return to step 3. If tol3 is greater than or equal to TOL3, the final wave impedance and seismic wavelet inversion results are obtained. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] Figure 1It is a flow chart of the present invention. DETAILED DESCRIPTION

[0037] The specific technical solution of the present invention is explained in conjunction with embodiments.

[0038] A blind inversion method based on construction and non-convex total variation joint constraints includes the following steps:

[0039] Step 1: Input seismic data and calculate the initial seismic wavelet and initial wave impedance model. The specific steps are as follows:

[0040] Step 1.1: Input seismic data and perform stacking processing on the seismic data to obtain post-stack seismic data Where n is the number of sampling points in each channel, and m is the number of seismic channels.

[0041] Step 1.2: Perform spectrum analysis on the seismic data, estimate the main frequency f0 of the seismic data, and calculate the zero-phase Ricker wavelet with the main frequency f0 as the initial seismic wavelet:

[0042]

[0043] Where a(t) is the amplitude of the Ricker wavelet at time t. Determine the duration of the wavelet t n and sampling interval nt, then the jth initial seismic wavelet j w0=[a(t0),a(t1),…,a(t tn / nt )].

[0044] Step 1.3: Use the initial seismic wavelet to convolve the reflection coefficient sequence calculated from the well logging curve to obtain a synthetic seismic record, compare the synthetic seismic record with the actual seismic record until the main stratigraphic positions match, and then make fine adjustments to maximize the correlation coefficient, thereby obtaining the well logging curve in the time domain;

[0045] Starting from the well logging curve, the stratigraphic layer is tracked along the seismic phase axis, and a two-dimensional profile is obtained by interpolating the one-dimensional well logging data along the layer. The interpolated logarithmic wave impedance is low-pass filtered with a cutoff frequency of 5 Hz to obtain the initial wave impedance model. right Taking the logarithm we get

[0046] Step 2: Calculate the seismic dip angle through post-stack seismic data and construct a rotation operator. The specific steps are as follows:

[0047] Step 2.1: Calculate the gradient structure tensor ST:

[0048]

[0049] Among them G ρrepresents the Gaussian kernel function, ρ is the standard deviation of the Gaussian kernel, and d represents the first-order derivative vector at the sampling point of the two-dimensional seismic data. x and d y are the first-order derivatives in the horizontal and vertical directions respectively, and the eigenvalue decomposition of ST is performed:

[0050]

[0051] Where λ1 and λ2 are eigenvalues, and v1 and v2 are the corresponding eigenvectors. The formula for calculating the local inclination is:

[0052]

[0053] where v 11 and v 12 are the two components of the eigenvector corresponding to the principal eigenvalue, and ε is a value close to zero.

[0054] Step 2.2: The rotation operator can be expressed as:

[0055]

[0056] D perp is the first-order difference operator perpendicular to the local structure direction, D parl is a first-order difference operator parallel to the local structure. cos and Q sin are diagonal matrices composed of the cosine and sine values ​​of the inclination angle θ:

[0057]

[0058] Dx and Dy are the first-order difference operators in the horizontal direction and the vertical direction respectively:

[0059]

[0060] Step 3: Establish the objective function of wave impedance inversion and solve it using the split-bregman algorithm. The specific steps include:

[0061] Step 3.1: Rearrange the seismic traces into 1-dimensional vectors in Two-dimensional seismic record The jth lane of

[0062] Rearrange the initial wave impedance model into a 1-dimensional vector in is the initial wave impedance model of the jth channel;

[0063] Step 3.2: Since this inversion algorithm is based on the self-excited and self-received seismic profile, we can invert the corresponding seismic wavelet for each channel and define the wavelet matrix of the jth channel: j W:

[0064]

[0065] When k = 0, each seismic wavelet is an initial wavelet, that is, j w= j w0.

[0066] Define the first-order difference operator:

[0067]

[0068] Define the single channel coefficient matrix j G= j WD; therefore, the multi-channel coefficient matrix G is expressed as:

[0069]

[0070] Step 3.3: Establish the objective function of wave impedance inversion:

[0071]

[0072] Where Z is the logarithmic impedance to be inverted, Φ α () is the smooth clipped absolute deviation (SCAD) penalty, R perp and R parl are auxiliary variables. μ and σ are regularization parameters.

[0073] R parl For example, the SCAD penalty operator Φ α (R parl ) is expressed as:

[0074]

[0075] Among them, r s is the vector R parl The sth element in , α is the SCAD penalty factor, and ζ is the shape tuning parameter.

[0076] Step 3.4: Convert the objective function into an unconstrained optimization problem:

[0077]

[0078] where b prep and b parl Ask about auxiliary variables separately.

[0079] Input parameters: μ, γ, σ, α and ξ, input iteration stop parameter tol1, initialize auxiliary variables Let k = 0, where the superscript k is the kth iteration, and use the split-bregman algorithm to solve the objective function, splitting the original problem into three sub-problems:

[0080]

[0081] Step 3.5: Solve the subproblems Update the logarithmic impedance Z, which is a linear inversion problem subject to strict quadratic regularization constraints. Therefore, the exact solution to the equation is:

[0082]

[0083] where Λ is the identity matrix.

[0084] Step 3.6: Solve the subproblem using the SCAD threshold operator Updated R perp for:

[0085]

[0086] Solving the subproblem using the SCAD threshold operator Updated R parl for:

[0087]

[0088] Step 3.7: Update auxiliary variable b perp and b parl :

[0089]

[0090] Step 3.8: Define the iteration stopping criterion TOL1:

[0091]

[0092] Determine the logarithmic impedance Z k+1 Whether the iteration stopping criteria are met: If tol1 is less than TOL1, set k = k + 1 and return to step 3.5. If tol1 is greater than or equal to TOL1, output Z k+1 .

[0093] Step 4: Establish the objective function of wavelet inversion and solve it using the split-bregman algorithm. The specific steps are:

[0094] Step 4.1: Define the coefficient matrix R of wavelet inversion, the reflection coefficient sequence of the jth channel Its Toblitz matrix formj R is:

[0095]

[0096] The multi-channel coefficient matrix is:

[0097]

[0098] Step 4.2: Arrange the initial seismic wavelet of each channel of the two-dimensional seismic profile into a one-dimensional vector

[0099] Step 4.3: Establish the objective function of wavelet inversion:

[0100]

[0101] in is the wavelet to be inverted, μ1, μ2, μ3, μ4 are regularization parameters, R1, R2 are auxiliary variables, and D1 is a first-order difference operator, which is formally equal to D y , D2 is the second-order difference operator, expressed as:

[0102]

[0103] Convert the constrained optimization problem into an unconstrained optimization problem:

[0104]

[0105] Among them, γ1 and γ2 are positive parameters, and b1 and b2 are auxiliary variables.

[0106] Given the regularization parameter, input the iteration stop parameter tol2, and initialize the auxiliary variables Split the original problem into 3 sub-problems:

[0107]

[0108] Step 4.4: Solve the subproblems Update multi-channel sub-wave vectors This is a linear inversion problem subject to strict quadratic regularization constraints. Therefore, the exact solution to the equation is:

[0109]

[0110] Step 4.5: Solve the subproblem using the soft threshold operator have to:

[0111]

[0112] Solve the subproblem using the soft threshold operator have to:

[0113]

[0114] Step 4.6: Update auxiliary variables b1, b2:

[0115]

[0116] Step 4.7: Define the iteration stopping criterion TOL2:

[0117]

[0118] Determine whether the sub-wave vector meets the iteration stopping criteria: If tol2 is less than TOL2, set k = k + 1 and return to step 4.4. If tol2 is greater than or equal to TOL2, output

[0119] Step 5: Determine whether the iteration criteria are met. If not, return to step 3. If so, output the final wave impedance and wavelet inversion results. The specific steps are:

[0120] Step 5.1: Define outer loop iteration criteria TOL3:

[0121]

[0122] Determine whether the external iteration stopping criteria are met: given tol3, if tol3 is less than TOL3, then let Return to step 3. If tol3 is greater than or equal to TOL3, then set

[0123] Among them, exp() is the exponential operation, and reshape() is to restore the vector to an n×m-dimensional matrix. is the inverted two-dimensional wave impedance profile, is the inverted multi-channel seismic wavelet matrix.

Claims

1. A blind inversion method based on construction and non-convex total variation joint constraints, characterized in that: include: Step 1: Input seismic data and calculate the initial seismic wavelet and initial wave impedance model; Step 2: Calculate the seismic dip angle through post-stack seismic data and construct a rotation operator; Step 3: Establish the objective function of wave impedance inversion, bring in the initial seismic wavelet, and use the split-bregman algorithm to solve and obtain the wave impedance model; Step 4: Establish the objective function of wavelet inversion, bring in the wave impedance model, and use the split-bregman algorithm to solve and obtain the seismic wavelet; Step 5: Determine whether the iteration criteria are met. If not, return to step 3. If so, output the final wave impedance and wavelet inversion results.

2. A blind inversion method based on construction and non-convex total variation joint constraints according to claim 1, characterized in that: The step 1 comprises the following: Input seismic data, and perform stacking processing on the seismic data to obtain post-stack seismic data; Perform spectrum analysis on the seismic data, estimate the main frequency of the seismic data, and calculate the zero-phase Ricker wavelet of the corresponding main frequency as the initial seismic wavelet: The initial seismic wavelet is convolved with the reflection coefficient sequence calculated from the well logging curve to obtain a synthetic seismic record, which is then compared with the actual seismic record until the main stratigraphic positions match, and then fine-tuned to maximize the correlation coefficient, thereby obtaining the well logging curve in the time domain; Starting from the logging curve, the stratigraphic layer is tracked along the seismic phase axis, and a two-dimensional profile is obtained by interpolating the one-dimensional logging data along the layer. The interpolated logarithmic wave impedance is low-pass filtered with a cutoff frequency of 5 Hz to obtain the initial wave impedance model, and its logarithm is calculated.

3. A blind inversion method based on construction and non-convex total variation joint constraints according to claim 1, characterized in that: The step 2 includes the following: The gradient structure tensor is calculated based on the seismic data, the local dip is calculated from the gradient structure tensor, and the rotation operator is constructed from the local dip: D perp is a first-order difference rotation operator perpendicular to the local structure direction, D parl is a first-order difference rotation operator parallel to the local structure. cos and Q sin are diagonal matrices consisting of the cosine and sine of the local inclination angle, respectively.

4. The blind inversion method based on construction and non-convex total variation joint constraints according to claim 1, characterized in that: The step 3 includes the following: Establish the objective function of wave impedance inversion: subjecttoR perp =D perp Z,R parl =D parl Z Where G = WD, W is the wavelet matrix, D is the 0.5 times first-order difference matrix, Z is the logarithmic impedance to be inverted, Φ α () is the non-convex smooth clipped absolute deviation (SCAD) penalty, R perp and R parl are auxiliary variables, μ and σ are regularization parameters; The split-bregman algorithm is used to solve the objective function and obtain the wave impedance model.

5. The blind inversion method based on construction and non-convex total variation joint constraints according to claim 1, characterized in that: The step 4 includes the following: Establish the objective function of wavelet inversion: in is the seismic sub-wave vector, R is the diagonal matrix composed of the Toeplitz matrix of each reflection coefficient sequence, is the initial seismic wavelet, D1 and D2 are the first-order difference matrix and the second-order difference matrix respectively. μ1, μ2, μ3, μ4 are regularization parameters; The split-bregman algorithm is used to solve the objective function and obtain the seismic wavelet.

6. A blind inversion method based on construction and non-convex total variation joint constraints according to claim 1, characterized in that: The step 5 includes the following: Define the iteration stopping criterion: The superscript k+1 indicates the result of the k+1th cycle in the wavelet inversion cycle; Given the iteration stop parameter tol3, determine whether the external iteration stop criterion is met: if tol3 is less than TOL3, then let Return to step 3. If tol3 is greater than or equal to TOL3, the final wave impedance and seismic wavelet inversion results are obtained.

Images