A blind inversion method based on construction and non-convex total variation joint constraint

Abstract

The application provides a blind inversion method based on a structure and non-convex total variation joint constraint, and belongs to the field of seismic inversion. The method comprises the following steps: inputting seismic data, determining a post-stack seismic profile to be inverted; 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; establishing an objective function of wave impedance inversion, and solving the objective function by using a split-bregman algorithm to output wave impedance data; establishing an objective function of wavelet inversion, and bringing the wave impedance data into the objective function, and solving the objective function by using the split-bregman algorithm to output a seismic wavelet; alternately updating the objective functions of the wave impedance and the wavelet, and outputting the inverted wave impedance and the seismic wavelet when an iteration criterion is met. The application provides a wave impedance and wavelet simultaneous inversion framework, and the introduced structure constraint and non-convex total variation regularization can effectively improve the precision of wave impedance and wavelet inversion.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of coal, oil and gas and other mineral resources seismic exploration, and particularly relates to a blind inversion method based on joint constraint of structure and non-convex total variation. BACKGROUND

[0002] Seismic exploration is an important method for coal, oil and gas resource exploration. Wave impedance inversion is the main means to obtain underground wave impedance information from seismic data. It converts seismic data to a higher resolution model space, so that an accurate wave impedance model can be constructed. The conventional wave impedance inversion method uses total variation (TV) regularization to invert wave impedance underground. The TV regularization constrains the horizontal and vertical derivatives of the model parameters using the L1 norm, thereby constructing the objective function. Zhang (2014) et al. used the L1 norm mismatch function, TV regularization and prior information constraint to establish a seismic inversion objective function through the Lagrange multiplier method, solving the inversion problem of abnormal values and discontinuous boundaries in seismic data. Gholami (2015) used TV regularization to recover a blocky wave impedance model from reflectivity data, improving the spatio-temporal correlation of model parameters. Wang et al. (2018) proposed that, compared with isotropic total variation (ITV), anisotropic total variation (ATV) better preserves the discontinuous interfaces between strata. Zhang et al. (2022) developed an impedance inversion method based on TV constraint guided by geological structure (GSGTV), improving resolution and spatial continuity.

[0003] However, the seismic wavelet is unknown, and in general, a theoretical Ricker wavelet or a statistical wavelet can be used for inversion, but this will bring some errors. Therefore, Gholami (2016) developed a blind inversion method aimed at simultaneously recovering wavelet and wave impedance information from seismic data. However, TV regularization cannot fully exploit sparsity and has inherent bias, which cannot accurately recover stratum boundary information. In addition, TV regularization does not consider geological structure and is not suitable for seismic inversion in complex areas. SUMMARY

[0004] The present application proposes a blind inversion method based on joint constraint of structure and non-convex total variation to address the shortcomings of existing blind inversion methods, mainly solving the problems of inherent bias in TV regularization, which cannot accurately recover stratum boundary information and is not suitable for complex areas.

[0005] To solve the above technical problems, the technical solution adopted by the present application is: a blind inversion method based on joint constraint of structure and non-convex total variation, comprising:

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

[0007] Step 2: Calculate the seismic dip angle using 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. The synthetic seismic record is compared with the actual seismic record until the main stratigraphic layers coincide with each other. Then, fine-tuning is performed to maximize the correlation coefficient, and then the well logging curve in the time domain is obtained.

[0015] Starting from the well logging curve, the stratigraphic layer is traced 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, 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 differential rotation operator perpendicular to the local structure direction, D parl It 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 tilt angle, respectively.

[0020] Step 3 includes the following:

[0021] Objective function of wave impedance inversion is established:

[0022]

[0023] subject to R perp = D perp Z, R parl = D parl Z

[0024] Wherein G = WD, W is wavelet matrix, D is 0.5 times first-order difference matrix, Z is the log impedance to be inverted, Φ α () is non-convex smooth clipping 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 the wave impedance model is obtained.

[0026] Step 4 includes the following:

[0027] Objective function of wave impedance inversion is established:

[0028]

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

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

[0031] Step 5 includes the following:

[0032] Define iteration stopping criterion:

[0033]

[0034] Wherein, the superscript k+1 represents the result of the k+1th cycle in the wavelet inversion cycle;

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

[0036] Figure 1Flowchart of the present invention. DETAILED DESCRIPTION

[0037] The specific technical solutions of the present invention are described with reference to the 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. Then make fine adjustments to maximize the correlation coefficient, and then obtain the time domain well logging curve.

[0045] Starting from the well logging curve, the formation 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 using 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 ρis the Gaussian kernel function, p is the standard deviation of the Gaussian kernel, and d represents the first derivative vector at the 2D seismic data sampling point. x and d y are the first derivatives in the horizontal and vertical directions, respectively, and ST is subjected to eigenvalue decomposition:

[0050]

[0051] where λ1 and λ2 are the eigenvalues, and v1 and v2 are the corresponding eigenvectors. The formula for calculating the local dip angle 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 represented as:

[0055]

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

[0057]

[0058] Dx and Dy are the first difference operators in the horizontal and vertical directions, 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 1D vectors where is the jth trace of the 2D seismic record ;

[0062] Rearrange the initial wave impedance model into 1D vectors where is the jth trace of the initial wave impedance model;

[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 a 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 smoothed 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 of 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] Solve 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 criterion is 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 and 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 2D 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, D1 is the 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 as follows:

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

[0121]

[0122] Determine whether the external iteration stopping criterion is 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 joint constraints of construction and non-convex total variation, characterized in that: include: Step 1: Input seismic data and calculate the initial seismic wavelet and initial wave impedance model; Step 2: Calculate the gradient structure tensor based on the seismic data, calculate the local dip angle from the gradient structure tensor, and then construct the rotation operator from the local dip angle: D x is the first-order differential rotation operator in the horizontal direction, D y is the first-order difference operator in the vertical direction, D perp is a first-order differential rotation operator perpendicular to the local structure direction, D parl is a first-order differential rotation operator parallel to the local structure, Q cos and Q sin are diagonal matrices consisting of the cosine and sine values ​​of the local tilt angle; Step 3: Substitute the initial seismic wavelet and 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 smoothed clipped absolute deviation (SCAD) penalty, Z0 is the low-frequency logarithmic impedance, 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. Step 4: Establish the objective function of wavelet inversion and bring the wave impedance model into it: in is the seismic wavelet vector, R is the diagonal matrix formed by 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; 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 includes 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. The synthetic seismic record is compared with the actual seismic record until the main stratigraphic layers coincide with each other. Then, fine-tuning is performed to maximize the correlation coefficient, and then the well logging curve in the time domain is obtained. Starting from the well logging curve, the stratigraphic layer is traced 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, and its logarithm is calculated.

3. The 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 represents the result of the k+1th cycle in the wavelet inversion cycle; Given the iteration stopping parameter tol3, determine whether the external iteration stopping 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