Wave impedance inversion method based on reweighted arc tangent total variation constraint

Abstract

The invention belongs to the technical field of seismic data processing and oil-gas exploration, provides a wave impedance inversion method based on a reweighted arc tangent total variation constraint, and aims to solve the problems of common pseudo-layer phenomenon, insufficient transverse continuity and easy falling into a local optimal solution in wave impedance inversion of a traditional L1 norm and Lp quasi-norm constraint method. According to the method, the reweighted arc tangent total variation constraint is innovatively introduced, the alternating direction multiplier method (ADMM) and the Sylvester equation are combined for solving, the weight matrix is dynamically adjusted, generation of pseudo layers is effectively reduced, and the transverse continuity and the overall precision of an inversion result are remarkably improved. In addition, the method has the advantages of high calculation efficiency, high stability and the like. In oil-gas exploration, the method can provide a more accurate and reliable wave impedance model, greatly enhances the recognition capability of oil-gas reservoirs, is suitable for the fields of seismic data processing and oil-gas resource development, and has important practical application value for improving the oil-gas exploration efficiency and accuracy.

Description

Technical Field

[0001] The invention relates to seismic wave impedance inversion, belongs to the technical field of oil and gas exploration and seismic inversion, and provides a wave impedance inversion method based on reweighted inverse tangent total variation constraint. Background Art

[0002] Seismic wave impedance inversion uses the mathematical and physical relationship between post-stack seismic records and wave impedance to establish a mathematical model, and uses the optimization method to solve the mathematical model to obtain seismic wave impedance from seismic records. It is an important means of oil and gas exploration. The seismic wave impedance inversion technology based on sparse constraints is an important method of seismic wave impedance inversion. It introduces sparse constraints as prior information into the forward model, successfully improving the accuracy of the inversion results.

[0003] The seismic impedance inversion method based on L1 norm constraint is an important method of seismic impedance inversion technology based on sparse constraints. This method uses L1 norm as a constraint term to construct a forward model, and iteratively solves the optimization problem of the forward model through the optimization method, and finally obtains the impedance inversion result. Li et al. first introduced the anisotropic total variation regularization based on L1 norm into the impedance inversion, and proved that this method can improve the inversion accuracy. He et al. used the weighted L1 norm constraint for single-channel impedance inversion and used the alternating direction multiplier method to solve it, proving the feasibility of this method in seismic impedance inversion. However, the L1 norm cannot fully exploit the sparsity and cannot reflect the relationship between channels, resulting in the inversion results having pseudo layers and poor lateral continuity, thereby reducing the accuracy. Wu et al. combined fractional differences with ATpV and used the fractional anisotropic total variation method to solve the problem that the inversion results caused by the L1 norm have pseudo layers and poor lateral continuity. However, the inversion method based on the Lp pseudo-norm will fall into local optimality during the solution process, resulting in poor recovery effect near pinch points and thin layers, thereby reducing the inversion accuracy. Therefore, an inversion method is needed to overcome the above problems. Summary of the invention

[0004] The purpose of the present invention is to solve the problems in the prior art of seismic wave impedance inversion method based on L1 norm and Lp pseudo-norm constraints, such as pseudo layers, poor lateral continuity, and easy falling into local optimum, so as to improve the accuracy of wave impedance inversion.

[0005] In order to achieve the above object, the present invention adopts the following technical means:

[0006] The present invention provides a wave impedance inversion method based on reweighted inverse tangent total variation constraint, comprising the following steps:

[0007] Step 1: Preprocess multi-channel seismic data and obtain an initial model profile of wave impedance, wherein the seismic data includes seismic records, well logging data and layer information;

[0008] Step 2: Calculate the initial logarithmic wave impedance according to the initial model section of the wave impedance, initialize the reweighted matrix, and establish a seismic wave impedance forward model based on the reweighted inverse tangent total variation constraint;

[0009] Step 3: Input the initial matrix of Lagrange multiplier terms and the initial matrix of its dual terms, and update the initial logarithmic wave impedance, the initial matrix of Lagrange multiplier terms and the initial matrix of its dual terms by combining the alternating direction multiplier method and the seismic wave impedance forward model constrained by Sylvester equation and inverse tangent total variation to obtain the updated logarithmic wave impedance;

[0010] Step 4: Determine whether the values ​​of the initial logarithmic wave impedance before and after the update meet the conditions. If so, update the updated logarithmic wave impedance and Lagrange multiplier and the dual term As the initial value, and go to step 3, otherwise, perform exponential operation based on the updated logarithmic wave impedance to obtain the inversion result of the wave impedance.

[0011] In the above scheme, step 1 comprises the following steps:

[0012] Step 1.1: Extract post-stack seismic records, well logging data and stratigraphic information from multi-channel seismic data, and perform low-pass filtering on the well logging data to obtain the initial model profile of the parameters to be inverted, that is, the initial model profile of the wave impedance.

[0013] In the above scheme, step 2 comprises the following steps:

[0014] Step 2.1: Based on the initial model section and logarithmic operation formula of the wave impedance, the initial logarithmic wave impedance L0 is obtained. The logarithmic operation formula is:

[0015] L=ln(Z) (1)

[0016] Wherein, Z represents wave impedance, Z0 is the initial wave impedance for solving the initial logarithmic wave impedance L0, L represents the logarithmic wave impedance matrix, and ln represents the natural logarithm operation;

[0017] Step 2.2: Use the longitudinal difference matrix D of the difference matrix D v , the transverse difference matrix D x Constructing the logarithmic wave impedance longitudinal differential L y and the logarithmic wave impedance transverse differential L x , and based on the logarithmic wave impedance longitudinal differential L y and the logarithmic wave impedance transverse differential L x Construct the reweighted inverse tangent total variation constraint term, the formula is:

[0018] L y =Dy L (2)

[0019] L x =LD x (3)

[0020] Longitudinal difference matrix D y The expression is as follows:

[0021]

[0022] Where m represents the number of sampling points of a single channel of logarithmic wave impedance;

[0023] Transverse difference matrix D x The expression is as follows:

[0024]

[0025] Where, t represents the number of channels of the logarithmic wave impedance profile;

[0026] The reweighted inverse tangent total variation constraint term is:

[0027]

[0028] Among them, a is the given inverse tangent total variation parameter, which determines the sparsity of the inverse tangent total variation, L x (i) represents the transverse differential impedance L x The number of the i-th channel, L y (j) represents the longitudinal differential impedance L y The number of the jth channel, M x With M y is the reweighting matrix, which controls the regularization coefficient of the horizontal and vertical differential terms update. The initial value is 1, and the update formula is as follows:

[0029]

[0030] ∈ is a constant used to prevent division by zero errors;

[0031] Step 2.3: Based on the given inverse tangent total variation parameter a, sparse constraint weight coefficient μ, initial model section weight coefficient ω and convolution matrix W, a seismic impedance forward model based on reweighted inverse tangent total variation constraint is constructed:

[0032]

[0033] Among them, ||*|| F represents the Frobenius norm, the matrix It refers to the set of all logarithmic wave impedances L that make the function J(L) reach its minimum value. WDL refers to the multiplication of the convolution matrix W, the difference matrix D, and the logarithmic wave impedance matrix L. ψ(*) means the reweighted inverse tangent total variation constraint term, and S0 means the actual seismic data.

[0034] The convolution matrix W is expressed as:

[0035]

[0036] Among them, w q represents the qth sampling point of the seismic wavelet, and m represents the number of sampling points of the logarithmic wave impedance single channel;

[0037] Step 2.4: Introduce the Lagrange multiplier term R based on the seismic impedance forward model x , R y and the dual term C x , C y , the objective function of seismic impedance inversion based on reweighted inverse tangent total variation constraint is obtained:

[0038]

[0039] Where a represents the total variation parameter of the inverse tangent, λ is the weight coefficient of the dual term, R x (i) represents the Lagrange multiplier term of the i-th logarithmic impedance profile, R y (j) represents the Lagrange multiplier of the jth sampling point of the logarithmic wave impedance single channel, and S0 represents the actual seismic record.

[0040] In the above scheme, step 3 comprises the following steps:

[0041] Step 3.1: Use the initial logarithmic wave impedance L0 as the initial model section L 1 , introduce the initial matrix of Lagrange multiplier term and its dual initial matrix in, represents the Lagrange multiplier term of the i-th logarithmic impedance profile after initialization, It represents the Lagrange multiplier term of the jth sampling point of the logarithmic wave impedance single channel after initialization, represents the dual term of the i-th logarithmic wave impedance profile after initialization, Represents the dual term of the jth sampling point of the logarithmic wave impedance single channel after initialization;

[0042] Step 3.2: Based on step 3.1, the alternating direction multiplier algorithm and the Sylvester equation solution method are used to solve the objective function to obtain the wave impedance logarithm L* update formula and update the wave impedance logarithm L*. The formula is as follows:

[0043]

[0044] Among them, vec represents the matrix vectorization operator, E n represents the n×n identity matrix, T represents the transpose of the matrix, represents the Kronecker product, and the matrices A, B, and C are expressed as:

[0045]

[0046] Among them, E m-1 It represents the unit matrix of size (m-1)×(m-1), and λ is the weight coefficient of the dual term.

[0047] In the above scheme, step 4 is to determine whether the values ​​of the initial logarithmic wave impedance before and after the update meet the condition || L * -L 1 ||2 / ||L * ||2>tol, where tol is the given error term value. If it is, the updated logarithmic wave impedance and Lagrange multiplier terms are and the dual term As the initial value, go to step 3, otherwise, perform exponential operation based on the updated logarithmic wave impedance to obtain the wave impedance inversion result Z r :

[0048] Z f =exp(L*) (16)

[0049] Here, exp represents the exponential operator.

[0050] In the above scheme, the Lagrange multiplier term and the dual term The steps to update are:

[0051] Calculate L 1 Logarithmic wave impedance transverse differential and log wave impedance longitudinal differential And update the Lagrange multiplier term according to the alternating direction multiplier method and the soft threshold method The calculation formula is as follows:

[0052]

[0053]

[0054] ∈ is a constant used to prevent division by zero errors, sign represents the sign function, and the expression is:

[0055]

[0056] Based on the wave impedance logarithm L* and the updated Lagrange multiplier term Update the dual term according to the alternating direction multiplier algorithm The calculation formula is as follows:

[0057]

[0058] Because the present invention adopts the above technical means, it has the following beneficial effects:

[0059] 1. The present invention adopts a reweighted inverse tangent total variation constraint, combined with the alternating direction multiplier method (ADMM) and the Sylvester equation for solution. The invention solves the technical problem in the prior art that the sparse constraints based on the L1 norm and the Lp pseudo-norm are prone to produce pseudo-layer phenomena in seismic wave impedance inversion, and the lateral continuity of the inversion results is poor, resulting in low inversion accuracy. By introducing the reweighted inverse tangent total variation constraint, the generation of pseudo-layers is effectively reduced, the lateral continuity of the inversion results is improved, and thus the accuracy of wave impedance inversion is significantly improved.

[0060] Second, by combining the alternating direction multiplier method (ADMM) with the solution of the Sylvester equation, the technical problem that the existing inversion methods are prone to fall into the local optimum during the iteration process, resulting in the inability to obtain the global optimal inversion results, is solved. The ADMM method can effectively decompose complex optimization problems, and combined with the solution of the Sylvester equation, it ensures that a solution closer to the global optimum can be found during the inversion process, avoiding the trap of the local optimum, thereby improving the reliability of the inversion results.

[0061] 3. Initialize the reweighted matrix and dynamically update the weights according to the difference of logarithmic wave impedance during the iteration process. This solves the technical problem that the weights are usually fixed and cannot be adjusted adaptively in the traditional sparse constraint method, resulting in inaccurate inversion results. Dynamically updating the reweighted matrix enables the constraint terms to be adaptively adjusted, better captures the sparse characteristics of wave impedance, and improves the precision and accuracy of the inversion.

[0062] Fourth, by introducing Lagrange multiplier terms and dual terms and updating these terms in combination with the soft threshold method, the problem of difficulty in optimizing constraint terms in the inversion process is solved, achieving a better balance between regularization and data fitting.

[0063] 5. Through the preprocessing of multi-channel seismic data and the construction of the initial model, the problem that the single-channel inversion method cannot fully utilize the information of multi-channel seismic data, resulting in poor lateral continuity and overall consistency of the inversion results, is solved. Through the joint processing of multi-channel seismic data, the correlation between each channel can be fully utilized, the lateral continuity and overall consistency of the inversion results are improved, thereby improving the overall quality of the inversion.

[0064] In summary, the present invention effectively solves the problems of pseudo-layer, poor lateral continuity and easy falling into local optimum in the wave impedance inversion in the prior art through technical means such as reweighted inverse tangent total variation constraint, alternating direction multiplier method, solving Sylvester equation, dynamically updating reweighted matrix, introducing Lagrange multiplier terms and dual terms, etc., significantly improves the accuracy and reliability of the inversion results, and has important application value. BRIEF DESCRIPTION OF THE DRAWINGS

[0065] Figure 1 It is a flowchart of the present invention;

[0066] Figure 2 is the seismic profile in the present invention, which is a post-stack seismic record;

[0067] Figure 3 It is a schematic diagram of the cross section of the initial model of multi-channel wave impedance in the present invention;

[0068] Figure 4 It is a schematic diagram of the multi-channel wave impedance inversion result in the present invention;

[0069] Figure 5 To extract the seismic single channel beside the measured well and analyze it as a single channel. DETAILED DESCRIPTION

[0070] The following is a detailed description of the embodiments of the present invention. Although the present invention will be described and illustrated in conjunction with some specific embodiments, it should be noted that the present invention is not limited to these embodiments. On the contrary, modifications or equivalent substitutions made to the present invention should all be included in the scope of the claims of the present invention.

[0071] In addition, in order to better illustrate the present invention, numerous specific details are given in the following specific embodiments. It will be understood by those skilled in the art that the present invention can also be implemented without these specific details.

[0072] Combine the following Figure 1-5 The present invention is described in detail.

[0073] like Figure 1 As shown, a wave impedance inversion method based on reweighted inverse tangent total variation constraint includes the following steps:

[0074] Step 1: Preprocess multi-channel seismic data and obtain an initial model profile of wave impedance, wherein the seismic data includes seismic records, well logging data and layer information;

[0075] Step 1.1: Extract post-stack seismic records, well logging data and stratigraphic information from multi-channel seismic data, and perform low-pass filtering on the well logging data to obtain the initial model section of the parameters to be inverted, that is, the initial model section of the wave impedance, such as Figure 3 shown.

[0076] Step 2: Calculate the initial logarithmic wave impedance according to the initial model section of the wave impedance, initialize the reweighted matrix, and establish a seismic wave impedance forward model based on the reweighted inverse tangent total variation constraint;

[0077] Step 2.1: Based on the initial model section and logarithmic operation formula of the wave impedance, the initial logarithmic wave impedance L0 is obtained. The logarithmic operation formula is:

[0078] L=ln(Z)(1)

[0079] Wherein, Z represents wave impedance, Z0 is the initial wave impedance for solving the initial logarithmic wave impedance L0, L represents logarithmic wave impedance, and ln represents natural logarithm operation;

[0080] Step 2.2: Use the longitudinal difference matrix D of the difference matrix D y , the transverse difference matrix D x Constructing the logarithmic wave impedance longitudinal differential L y and the logarithmic wave impedance transverse differential L x , and based on the logarithmic wave impedance longitudinal differential L y and the logarithmic wave impedance transverse differential L x Construct the reweighted inverse tangent total variation constraint term, the formula is:

[0081] L y =D y L (2)

[0082] L x =LD x (3)

[0083] Longitudinal difference matrix D y The expression is as follows:

[0084]

[0085] Where m represents the number of sampling points of a single channel of logarithmic wave impedance;

[0086] Transverse difference matrix D x The expression is as follows:

[0087]

[0088] Where, t represents the number of channels of the logarithmic wave impedance profile;

[0089] The reweighted inverse tangent total variation constraint term is:

[0090]

[0091] Among them, atan represents the inverse tangent function, a is the given inverse tangent total variation parameter, which determines the sparsity of the inverse tangent total variation, and L x (i) represents the transverse differential impedance L x The number of the i-th channel, L y (j) represents the longitudinal differential impedance L y The number of the jth channel, M x With M y is the reweighting matrix, which controls the regularization coefficient of the horizontal and vertical differential terms update. The initial value is 1, and the update formula is as follows:

[0092]

[0093] ∈ is a very small constant, the value of ∈ is 0.00045, which is used to prevent division by zero errors.

[0094] Step 2.3: Based on the given inverse tangent total variation parameter a, sparse constraint weight coefficient μ, initial model section weight coefficient ω and convolution matrix W, a seismic impedance forward model based on reweighted inverse tangent total variation constraint is constructed:

[0095]

[0096] Among them, ||*|| F represents the Frobenius norm, the matrix It refers to the set of all logarithmic wave impedances L that make the function J(L) reach its minimum value, and WDL refers to the matrix multiplication of the convolution matrix W, the difference matrix D, and the logarithmic wave impedance matrix L;

[0097] The convolution matrix W is expressed as:

[0098]

[0099] Among them, w q represents the qth sampling point of the seismic wavelet;

[0100] Step 2.4: Introduce the Lagrange multiplier term R based on the seismic impedance forward model x , R y and the dual term C x , C y , the objective function of seismic impedance inversion based on reweighted inverse tangent total variation constraint is obtained:

[0101]

[0102] Among them, λ is the dual term weight coefficient, R x (i) represents the Lagrange multiplier term of the i-th logarithmic impedance profile, R y(j) represents the Lagrange multiplier term of the jth sampling point of the logarithmic wave impedance single channel.

[0103] Step 3: Input the initial matrix of Lagrange multiplier terms and the initial matrix of its dual terms, and update the initial logarithmic wave impedance, the initial matrix of Lagrange multiplier terms and the initial matrix of its dual terms by combining the alternating direction multiplier method and the seismic wave impedance forward model with the Svlvester equation and the inverse tangent total variation constraint to obtain the updated logarithmic wave impedance;

[0104] The specific steps are:

[0105] Step 3.1: Use the initial logarithmic wave impedance L0 as the initial model section L 1 , introduce the initial matrix of Lagrange multiplier term and its dual initial matrix in, represents the Lagrange multiplier term of the i-th logarithmic impedance profile after initialization, It represents the Lagrange multiplier term of the jth sampling point of the logarithmic wave impedance single channel after initialization, represents the dual term of the i-th logarithmic wave impedance profile after initialization, Represents the dual term of the jth sampling point of the logarithmic wave impedance single channel after initialization;

[0106] Step 3.2: Based on step 3.1, the alternating direction multiplier algorithm and the Sylvester equation solution method are used to solve the objective function to obtain the wave impedance logarithm L* update formula and update the wave impedance logarithm L*. The formula is as follows:

[0107]

[0108] Among them, vec represents the matrix vectorization operator, E n represents the n×n identity matrix, T represents the transpose of the matrix, represents the Kronecker product, and the matrices A, B, and C are expressed as:

[0109]

[0110] Among them, E m-1 Represents the identity matrix of size (m-1)×(m-1).

[0111] Step 4: Determine whether the values ​​of the initial logarithmic wave impedance before and after the update meet the conditions. If so, update the updated logarithmic wave impedance and Lagrange multiplier and the dual term As the initial value, go to step 3, otherwise, perform exponential operation based on the updated logarithmic wave impedance to obtain the inversion result of the wave impedance, such as Figure 4 shown.

[0112] Step 4 is to determine whether the values ​​of the initial logarithmic wave impedance before and after the update meet the condition || L * -L 1 ||2 / ||L * ||2>tol, where tol is the given error term value. If it is, the updated logarithmic wave impedance and Lagrange multiplier terms are and the dual term As the initial value, go to step 3, otherwise, perform exponential operation based on the updated logarithmic wave impedance to obtain the wave impedance inversion result Z r :

[0113] Z r =exp(L * ) (16)

[0114] Here, exp represents the exponential operator.

[0115] Lagrange multiplier term and the dual term The steps to update are:

[0116] Calculate L 1 Logarithmic wave impedance transverse differential and log wave impedance longitudinal differential And update the Lagrange multiplier term according to the alternating direction multiplier method and the soft threshold method The calculation formula is as follows:

[0117]

[0118] ∈ is a small constant used to prevent division by zero errors. Sign represents the sign function. The expression is:

[0119]

[0120] Based on the wave impedance logarithm L* and the updated Lagrange multiplier term Update the dual term according to the alternating direction multiplier algorithm The calculation formula is as follows:

[0121]

[0122] The present invention also provides a seismic wave impedance inversion system based on logarithmic total variation sparse constraints, comprising:

[0123] Preprocessing module: preprocesses multi-channel seismic data and obtains the initial model profile of wave impedance, where the seismic data includes seismic records, well logging data and layer information;

[0124] Seismic wave impedance forward model building module: calculate the initial logarithmic wave impedance according to the initial model section of wave impedance, initialize the reweighted matrix, and establish the seismic wave impedance forward model based on the reweighted inverse tangent total variation constraint;

[0125] Update module: input the initial matrix of Lagrange multiplier terms and the initial matrix of its dual terms, and update the initial logarithmic wave impedance, the initial matrix of Lagrange multiplier terms and the initial matrix of its dual terms by combining the alternating direction multiplier method and the seismic wave impedance forward model with Sylvester equation and reweighted inverse tangent total variation constraint to obtain the updated logarithmic wave impedance;

[0126] Judgment module: judge whether the values ​​of the initial logarithmic wave impedance before and after the update meet the conditions. If so, the updated logarithmic wave impedance and Lagrange multiplier items are and the dual term As the initial value, and go to step 3, otherwise, perform exponential operation based on the updated logarithmic wave impedance to obtain the inversion result of the wave impedance.

[0127] Furthermore, the specific implementation steps of the preprocessing module are:

[0128] Step 1.1: Extract post-stack seismic records, well logging data and stratigraphic information from multi-channel seismic data, and perform low-pass filtering on the well logging data to obtain the initial model section of the parameters to be inverted, that is, the initial model section of the wave impedance, such as Figure 3 shown.

[0129] Furthermore, the specific implementation steps of the seismic wave impedance forward model establishment module are:

[0130] Step 2.1: Based on the initial model section and logarithmic operation formula of the wave impedance, the initial logarithmic wave impedance L0 is obtained. The logarithmic operation formula is:

[0131] L=ln(Z) (1)

[0132] Wherein, Z represents wave impedance, Z0 is the initial wave impedance for solving the initial logarithmic wave impedance L0, L represents logarithmic wave impedance, and ln represents natural logarithm operation;

[0133] Step 2.2: Use the longitudinal difference matrix D y , the transverse difference matrix D x Constructing the logarithmic wave impedance longitudinal differential L y and the logarithmic wave impedance transverse differential L x , and based on the logarithmic wave impedance longitudinal differential L y and the logarithmic wave impedance transverse differential L x Construct the reweighted inverse tangent total variation constraint term, the formula is:

[0134] Ly =D y L (2)

[0135] L x =LD x (3)

[0136] Longitudinal difference matrix D y The expression is as follows:

[0137]

[0138] Where m represents the number of sampling points of a single channel of logarithmic wave impedance;

[0139] Transverse difference matrix D x The expression is as follows:

[0140]

[0141] Where, t represents the number of channels of the logarithmic wave impedance profile;

[0142] The reweighted inverse tangent total variation constraint term is:

[0143]

[0144] Among them, atan represents the inverse tangent function, a is the given inverse tangent total variation parameter, which determines the sparsity of the inverse tangent total variation, and L x (i) represents the transverse differential impedance L x The number of the i-th channel, L y (j) represents the longitudinal differential impedance L y The number of the jth channel, M x With M y is the reweighting matrix, which controls the regularization coefficient of the horizontal and vertical differential terms update. The initial value is 1, and the update formula is as follows:

[0145]

[0146] ∈ is a small constant used to prevent division by zero errors.

[0147] Step 2.3: Based on the given inverse tangent total variation parameter a, sparse constraint weight coefficient μ, initial model section weight coefficient α and convolution matrix W, a seismic impedance forward model based on reweighted inverse tangent total variation constraint is constructed:

[0148]

[0149] Among them, ||*|| F represents the Frobenius norm, the matrix It refers to the set of all logarithmic wave impedances L that make the function J(L) reach its minimum value, and WDL refers to the multiplication of the convolution matrix W, the difference matrix D, and the logarithmic wave impedance matrix L;

[0150] The convolution matrix W is expressed as:

[0151]

[0152] Among them, w q represents the qth sampling point of the seismic wavelet;

[0153] Step 2.4: Introduce the Lagrange multiplier term R based on the seismic impedance forward model x , R y and the dual term C x , C y , the objective function of seismic impedance inversion based on reweighted inverse tangent total variation constraint is obtained:

[0154]

[0155] Among them, λ is the dual term weight coefficient, R x (i) represents the Lagrange multiplier term of the i-th logarithmic impedance profile, R y (j) represents the Lagrange multiplier term of the jth sampling point of the logarithmic wave impedance single channel.

[0156] Furthermore, the specific implementation steps of the update module are:

[0157] Step 3.1: Use the initial logarithmic wave impedance L0 as the initial model section L 1 , introduce the initial matrix of Lagrange multiplier term and its dual initial matrix in, represents the Lagrange multiplier term of the i-th logarithmic impedance profile after initialization, It represents the Lagrange multiplier term of the jth sampling point of the logarithmic wave impedance single channel after initialization, represents the dual term of the i-th logarithmic wave impedance profile after initialization, Represents the dual term of the jth sampling point of the logarithmic wave impedance single channel after initialization;

[0158] Step 3.2: Based on step 3.1, the alternating direction multiplier algorithm and the Sylvester equation solution method are used to solve the objective function to obtain the wave impedance logarithm L* update formula and update the wave impedance logarithm L*. The formula is as follows:

[0159]

[0160] Among them, vec represents the matrix vectorization operator, En represents the n×n identity matrix, T represents the transpose of the matrix, represents the Kronecker product, and the matrices A, B, and C are expressed as:

[0161]

[0162] Among them, E m-1 Represents the identity matrix of size (m-1)×(m-1).

[0163] Judgment module: judge whether the values ​​of the initial logarithmic wave impedance before and after the update meet the conditions. If so, the updated logarithmic wave impedance and Lagrange multiplier items are and the dual term As the initial value, go to step 3, otherwise, perform exponential operation based on the updated logarithmic wave impedance to obtain the inversion result of the wave impedance, such as Figure 4 shown.

[0164] Determine whether the initial logarithmic wave impedance value before and after update meets the condition || L * -L 1 ||2 / ||L * ||2>tol, where tol is the given error term value. If it is, the updated logarithmic wave impedance and Lagrange multiplier terms are and the dual term As the initial value, go to step 3, otherwise, perform exponential operation based on the updated logarithmic wave impedance to obtain the wave impedance inversion result Z r :

[0165] Z r =exp(L*) (16)

[0166] Here, exp represents the exponential operator.

[0167] Lagrange multiplier term and the dual term The steps to update are:

[0168] Calculate L 1 Logarithmic wave impedance transverse differential and log wave impedance longitudinal differential And update the Lagrange multiplier term according to the alternating direction multiplier method and the soft threshold method The calculation formula is as follows:

[0169]

[0170] ∈ is a small constant used to prevent division by zero errors. Sign represents the sign function. The expression is:

[0171]

[0172] Based on the wave impedance logarithm L* and the updated Lagrange multiplier term Update the dual term according to the alternating direction multiplier algorithm The calculation formula is as follows:

[0173]

[0174] In summary, a seismic impedance inversion method based on reweighted inverse tangent total variation constraint is used to invert the impedance of multi-channel seismic profiles (two-dimensional seismic data obtained from seismic data), such as Figure 2-5 As shown in the figure, the horizontal axis is the number of seismic traces (Distance / Trace), and the vertical axis is the number of sampling points (Depth / Sampling point). Compared with the initial model, the inversion result of the wave impedance profile can better conform to the trend of the seismic record, which proves the correctness of this method. Figure 4 The inversion results have very good continuity, such as Figure 5 As shown, True Model represents the theoretical model, Initial Model represents the initial model, and Inversion Results represent the inversion results. The inversion results are very close to the real wave impedance and do not cause the pseudo layer phenomenon, which further proves the correctness of the method. The present invention replaces the L1 norm in the L1 norm sparse constraint with a reweighted inverse tangent total variation sparse constraint, and proposes a wave impedance inversion method based on the reweighted inverse tangent total variation constraint in combination with the alternating direction multiplier method. The initial model is inverted and updated using the inversion method mentioned above, and the optimal inversion result is output through repeated iterations. The present invention solves the problems in the prior art of seismic wave impedance inversion methods based on L1 norm and Lp pseudo-norm constraints, such as pseudo layers in the inversion results, poor lateral continuity, and easy to fall into local optimality, thereby improving the accuracy of wave impedance inversion.

[0175] The above are only representative embodiments of the present invention in many specific application scopes, and do not constitute any limitation on the protection scope of the present invention. Any technical solutions formed by transformation or equivalent replacement fall within the protection scope of the present invention.

Claims

1. A method for inversion of wave impedance based on reweighted inverse tangent total variation constraint, characterized in that: The steps include: Step 1: Preprocess multi-channel seismic data and obtain an initial model profile of wave impedance, wherein the seismic data includes seismic records, well logging data and layer information; Step 2: Calculate the initial logarithmic wave impedance according to the initial model section of the wave impedance, initialize the reweighted matrix, and establish a seismic wave impedance forward model based on the reweighted inverse tangent total variation constraint; Step 3: Input the initial matrix of Lagrange multiplier terms and the initial matrix of its dual terms, and update the initial logarithmic wave impedance, the initial matrix of Lagrange multiplier terms and the initial matrix of its dual terms by combining the alternating direction multiplier method and the seismic wave impedance forward model constrained by Sylvester equation and inverse tangent total variation to obtain the updated logarithmic wave impedance; Step 4: Determine whether the values ​​of the initial logarithmic wave impedance before and after the update meet the conditions. If so, update the updated logarithmic wave impedance and Lagrange multiplier and the dual term As the initial value, go to step 3, otherwise, perform exponential operation based on the updated logarithmic wave impedance to obtain the inversion result of the wave impedance.

2. The method of wave impedance inversion based on reweighted inverse tangent total variation constraint according to claim 1, characterized in that: The step 1 comprises the following steps: Step 1.1: Extract post-stack seismic records, well logging data and stratigraphic information from multi-channel seismic data, and perform low-pass filtering on the well logging data to obtain the initial model profile of the parameters to be inverted, that is, the initial model profile of the wave impedance.

3. The method of wave impedance inversion based on reweighted inverse tangent total variation constraint according to claim 1, characterized in that: The step 2 comprises the following steps: Step 2.1: Based on the initial model section and logarithmic operation formula of the wave impedance, the initial logarithmic wave impedance L0 is obtained. The logarithmic operation formula is: L=ln(Z)(1) Wherein, Z represents wave impedance, Z0 is the initial wave impedance for solving the initial logarithmic wave impedance L0, L represents the logarithmic wave impedance matrix, and ln represents the natural logarithm operation; Step 2.2: Use the longitudinal difference matrix D of the difference matrix D y , the transverse difference matrix D x Constructing the logarithmic wave impedance longitudinal differential L y and the logarithmic wave impedance transverse differential L x , and based on the logarithmic wave impedance longitudinal differential L y and the logarithmic wave impedance transverse differential L x Construct the reweighted inverse tangent total variation constraint term, the formula is: L y =D y L (2) L x =LD x (3) Longitudinal difference matrix D y The expression is as follows: Where m represents the number of sampling points of a single channel of logarithmic wave impedance; Transverse difference matrix D x The expression is as follows: Where, t represents the number of channels of the logarithmic wave impedance profile; The reweighted inverse tangent total variation constraint term is: Among them, a is the given inverse tangent total variation parameter, which determines the sparsity of the inverse tangent total variation, L x (i) represents the transverse differential impedance L x The number of the i-th channel, L y (j) represents the longitudinal differential impedance L y The number of the jth channel, M x With M y is the reweighting matrix, which controls the regularization coefficient of the horizontal and vertical differential terms update. The initial value is 1, and the update formula is as follows: ∈ is a constant used to prevent division by zero errors; Step 2.3: Based on the given inverse tangent total variation parameter a, sparse constraint weight coefficient μ, initial model section weight coefficient ω and convolution matrix W, a seismic impedance forward model based on reweighted inverse tangent total variation constraint is constructed: Among them, ||*|| F represents the Frobenius norm, the matrix It refers to the set of all logarithmic wave impedances L that make the function J(L) reach its minimum value. WDL refers to the multiplication of the convolution matrix W, the difference matrix D, and the logarithmic wave impedance matrix L. ψ(*) means the reweighted inverse tangent total variation constraint term, and S0 means the actual seismic data. The convolution matrix W is expressed as: Among them, w q represents the qth sampling point of the seismic wavelet, and m represents the number of sampling points of the logarithmic wave impedance single channel; Step 2.4: Introduce the Lagrange multiplier term R based on the seismic impedance forward model x , R y and the dual term C x , C y , the objective function of seismic impedance inversion based on reweighted inverse tangent total variation constraint is obtained: Where a represents the total variation parameter of the inverse tangent, λ is the weight coefficient of the dual term, R x (i) represents the Lagrange multiplier term of the i-th logarithmic impedance profile, R y (j) represents the Lagrange multiplier of the jth sampling point of the logarithmic wave impedance single channel, and S0 represents the actual seismic record.

4. The method of wave impedance inversion based on reweighted inverse tangent total variation constraint according to claim 3, characterized in that: The step 3 comprises the following steps: Step 3.1: Use the initial logarithmic wave impedance L0 as the initial model section L 1 , introduce the initial matrix of Lagrange multiplier term and its dual initial matrix in, represents the Lagrange multiplier term of the i-th logarithmic impedance profile after initialization, It represents the Lagrange multiplier term of the jth sampling point of the logarithmic wave impedance single channel after initialization, represents the dual term of the i-th logarithmic wave impedance profile after initialization, Represents the dual term of the jth sampling point of the logarithmic wave impedance single channel after initialization; Step 3.2: Based on step 3.1, the alternating direction multiplier algorithm and the Sylvester equation solution method are used to solve the objective function to obtain the wave impedance logarithm L* update formula and update the wave impedance logarithm L*. The formula is as follows: Among them, vec represents the matrix vectorization operator, E n represents the n×n identity matrix, T represents the transpose of the matrix, represents the Kronecker product, and the matrices A, B, and C are expressed as: Among them, E m-1 It represents the unit matrix of size (m-1)×(m-1), and λ is the weight coefficient of the dual term.

5. The method of wave impedance inversion based on reweighted inverse tangent total variation constraint according to claim 4, characterized in that: Step 4 is to determine whether the values ​​of the initial logarithmic wave impedance before and after the update meet the condition || L * -L 1 ||2 / ||L * ||2>tol, where tol is the given error term value. If it is, the updated logarithmic wave impedance and Lagrange multiplier terms are and the dual term As the initial value, go to step 3, otherwise, perform exponential operation based on the updated logarithmic wave impedance to obtain the wave impedance inversion result Z r : Z r =exp(L * ) (16) Here, exp represents the exponential operator.

6. The method of wave impedance inversion based on reweighted inverse tangent total variation constraint according to claim 5, characterized in that: Lagrange multiplier term and the dual term The steps to update are: Calculate L 1 Logarithmic wave impedance transverse differential and log wave impedance longitudinal differential And update the Lagrange multiplier term according to the alternating direction multiplier method and the soft threshold method The calculation formula is as follows: ∈ is a constant used to prevent division by zero errors, sign represents the sign function, and the expression is: Based on the wave impedance logarithm L* and the updated Lagrange multiplier term Update the dual term according to the alternating direction multiplier algorithm The calculation formula is as follows:

Images