Seismic wave impedance inversion method and system with minimum maximum concave high-order total variation constraint

Abstract

This invention discloses a seismic impedance inversion method and system based on minimum-maximum concavity higher-order total variation constraints. The method includes: preprocessing seismic and well logging data to construct an inversion objective function incorporating minimum-maximum concavity higher-order total variation constraints; and iteratively updating the natural logarithm of the impedance using the alternating direction multiplier method, thereby achieving high-precision seismic impedance inversion. The minimum-maximum concavity is used to describe high-frequency sparse information, and a higher-order total variation is employed to describe lateral and longitudinal differential information, effectively improving the inversion accuracy in structurally complex regions and suppressing artifacts in boundary regions. The advantages of this invention are: significantly improved resolution in regions with complex impedance structures; low error in the inversion results; and further enhanced resolution and accuracy of seismic impedance inversion.

Description

Technical Field

[0001] This invention relates to the fields of seismic inversion and oil and gas reservoir prediction, and in particular to a seismic impedance inversion method and system using minimum maximum concavity higher-order total variation constraints. Background Technology

[0002] Seismic impedance inversion is an important method for predicting hydrocarbon reservoirs. It establishes forward equations and an objective function based on the mathematical and physical relationship between post-stack seismic records and seismic impedance, and then uses optimization methods to solve the objective function to obtain the optimal solution for the seismic impedance. Sparse constraint-based seismic impedance inversion is a crucial method in seismic inversion, as it improves the accuracy of inversion results by introducing sparse constraints into the forward equations.

[0003] The traditional total variation constraint-based seismic impedance inversion method is an important approach to seismic impedance inversion techniques based on sparsity constraints. This method constructs a forward model and objective function using first-order total variation constraint terms, and iteratively solves the problem using the alternating direction multiplier method to obtain the impedance inversion result. This method can improve inversion accuracy. However, its ability to describe sparsity is weak, leading to low inversion accuracy in structurally complex regions and artifacts in boundary regions, thus reducing overall accuracy. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention provides a method and system for seismic impedance inversion with minimum maximum concavity higher-order total variational constraints.

[0005] To achieve the above-mentioned objectives, the technical solution adopted by the present invention is as follows:

[0006] A seismic impedance inversion method with minimum-maximum concavity higher-order total variational constraints includes the following steps:

[0007] S1. Preprocess the seismic data to obtain the initial data for seismic inversion, including: finite angle or post-stack seismic data, impedance data corresponding to the target stratigraphic segment, seismic wavelet data of the target layer, and the initial model for seismic inversion;

[0008] S2. Construct the seismic impedance inversion objective function constrained by minimum and maximum concavity higher-order total variational constraints, including: constructing a linear relationship model between reflection coefficient and natural logarithm of wave impedance based on data fidelity terms; constructing sparse constraint terms based on minimum and maximum concavity higher-order total variational constraints; and constructing the seismic impedance inversion objective function based on initial model constraint terms.

[0009] S3. Iteratively update the wave impedance logarithm using the alternating direction multiplier method, calculate and update the natural logarithm of the wave impedance, and update the Lagrange multipliers and dual terms.

[0010] S4. Calculate the inversion result of seismic impedance according to the conditions, and determine whether the values ​​of the natural logarithm of the wave impedance before and after the update meet the iterative update conditions. If so, take the updated natural logarithm of the wave impedance, Lagrange multipliers and dual terms as initial values ​​and transfer to the inversion iterative update module. If not, perform exponential operation on the updated natural logarithm of the wave impedance to obtain the inversion wave impedance result.

[0011] S5. Output and store the final inverted wave impedance results for easy subsequent analysis and application.

[0012] Further, step S1 includes the following sub-steps:

[0013] S101. Stack pre-stack seismic data according to the incident angle to obtain seismic data with limited angles or post-stack seismic data, referred to as seismic data, and the seismic data is denoted as S;

[0014] S102. Obtain the impedance data corresponding to the target formation section using well logging data. The impedance at a single well logging point is denoted as Z0.

[0015] S103. Using the seismic data from the wellbore and the target layer impedance at the well logging point, obtain the target layer seismic wavelet data w, w = [w1, w2, ..., w p ], where the length of the seismic wavelet data is p, w i (i = 1, 2, ..., p) represents the wavelet data value at the i-th sampling time;

[0016] S104. Based on the target formation stratigraphic data, perform low-pass filtering on the well logging data and extrapolate the values ​​according to the stratigraphic position to obtain the initial model profile of the wave impedance of the target stratigraphic segment.

[0017] Further, step S2 includes the following sub-steps:

[0018] S201. Calculate the initial natural logarithm of the wave impedance L0. The formula for calculating the natural logarithm is:

[0019] L0 = ln(Z0)

[0020] Where Z0 represents the initial wave impedance, L0 represents the natural logarithmic wave impedance, and ln represents the natural logarithmic operation formula;

[0021] S202. Construct a linear calculation model for the reflection coefficient R and its natural logarithm of wave impedance L, the calculation formula is as follows:

[0022] R = DL

[0023] In the above formula, R is the reflection coefficient, L is the natural logarithm of the wave impedance, and D is the difference matrix.

[0024] L = ln(Z)

[0025]

[0026] m represents the total number of sampling points for the logarithm L of the wave impedance;

[0027] S203. Construct the minimum and maximum concavity higher-order total variation constraints. The minimum and maximum concavity higher-order total variation constraints are calculated as follows:

[0028]

[0029] Among them, D xx D represents the second-order difference in the transverse direction. yy Let D represent the longitudinal second-order difference, γ represent the minimum-maximum concavity smoothing constraint parameter (whose value needs to be estimated in advance), ν represent the Moro envelope function factor, ||·||2 represents the L2 norm calculation formula, and ||·||1 represents the L1 norm calculation formula, where D xx L and LD yy The calculation formulas are as follows:

[0030] (D xx L) s,t =2L s,t -L s-1,t -L s+1,t

[0031] (LD yy ) s,t =2L s,t -L s,t-1 -L s,t+1

[0032] S204. Construct a seismic impedance forward model using minimum-maximum concavity higher-order total variational constraints:

[0033]

[0034] In the above formula, S represents seismic data, W represents the seismic wavelet convolution matrix, which is composed of the seismic wavelet data w described in step S103. The calculation method of the minimum-maximum concavity higher-order total variation constraint term is the same as in step S203. λ represents the sparse constraint weight coefficient, the value of which needs to be specified in advance. α represents the initial model constraint weight coefficient, the value of which needs to be specified in advance. The structure of the wavelet convolution matrix W is as follows:

[0035]

[0036] p represents the length of the seismic wavelet w.

[0037] S205. Construct an augmented Lagrangian function in the seismic wave impedance forward model to transform the constrained optimization problem into an unconstrained optimization problem:

[0038]

[0039] In the above formula, R xx and R yy Represents the Lagrange multiplier, C xx and C yy For the dual term, β is a parameter greater than 0.

[0040] Further, step S3 includes the following sub-steps:

[0041] S301. Use the logarithm of the initial model wave impedance L0 as the natural logarithm of the wave impedance to iteratively calculate the initial value L. 0 Take the initial value R of the Lagrange multipliers. xx =R yy =0 and the initial value of the dual term C xx =C yy =0;

[0042] S302. Update the natural logarithm L of the wave impedance using the solution to the Sylvester equation. k+1 The updated calculation formula is as follows:

[0043]

[0044] Among them, L k+1 Let represent the (k+1)th iteration of the logarithm of wave impedance, vec denotes the transformation of the matrix into a column vector, E represents the identity matrix, and T represents the transpose of the matrix. This represents the Kronecker product operation, where A, B, and C are represented as follows:

[0045]

[0046] S303. Update the Lagrange multiplier R xx and R yy The updated calculation formula is as follows:

[0047]

[0048] in, and This represents the result of the Lagrange multipliers in the (k+1)th iteration, where sign represents the sign function. and The calculation formula is:

[0049]

[0050] Among them, L k+1 This represents the (k+1)th iteration of the logarithm of the wave impedance. and This represents the result of the Lagrange multipliers in the k-th iteration. and This represents the result of the dual term in the k-th iteration;

[0051] S304. Update dual term C xx and C yy The updated calculation formula is as follows:

[0052]

[0053] in and This represents the result of the dual term in the (k+1)th iteration. and L represents the result of the dual term in the k-th iteration. k+1 This represents the (k+1)th iteration of the logarithm of the wave impedance. and This represents the result of the Lagrange multiplier subterms in the (k+1)th iteration.

[0054] Further, step S4 includes the following sub-steps:

[0055] S401. Update the logarithm L of the seismic impedance in the (k+1)th iteration. k+1 ;

[0056] S402. If ||L k+1 -L k ||2 / ||L k ||2>ε, where ε=5×10 -5 If the condition is met, return to step S302 for a loop; otherwise, calculate the wave impedance Z. The formula for calculating wave impedance Z is:

[0057] Z = exp(L k+1 )

[0058] Here, exp represents the natural exponent operation.

[0059] This invention also discloses a seismic impedance inversion system, which can be used to implement the above-mentioned seismic impedance inversion method, specifically including:

[0060] Data preprocessing module: Used to preprocess seismic data, including acquiring post-stack seismic record data, well logging data, seismic wavelets, and initial seismic inversion models.

[0061] Inversion objective function construction module: Creates a seismic wave impedance inversion calculation objective function based on data fidelity terms, minimum and maximum concavity high-order total variation constraints, and initial model constraints.

[0062] Inversion Iterative Update Module: The natural logarithm of wave impedance is updated based on the alternating direction multiplier method and the objective function of seismic wave impedance inversion, resulting in the updated natural logarithm of wave impedance.

[0063] Inversion Result Judgment Module: Determines whether the values ​​of the natural logarithm of wave impedance before and after the update meet the iterative update conditions. If so, the updated natural logarithm of wave impedance, Lagrange multipliers, and dual terms are used as initial values, and the module is transferred to the inversion iterative update module. If not, the updated natural logarithm of wave impedance is subjected to exponential operation to obtain the inversion wave impedance result.

[0064] Results Output Module: Outputs and stores the final inverted wave impedance results, including generating visualizations and detailed reports, ensuring that all output data and charts meet the expected format and accuracy requirements.

[0065] The present invention also discloses a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the above-described seismic wave impedance inversion method.

[0066] The present invention also discloses a computer-readable storage medium storing a computer program thereon, which, when executed by a processor, implements the above-described seismic wave impedance inversion method.

[0067] Compared with the prior art, the advantages of the present invention are as follows:

[0068] 1. By adopting the minimum-maximum concavity high-order total variation constraint method, the minimum-maximum concavity is used to describe high-frequency sparse information, which can more accurately capture subtle features and changes in seismic data and improve the accuracy of inversion.

[0069] 2. Using higher-order total variation to describe the lateral and longitudinal difference information can better preserve and reflect the spatial variation patterns of seismic data, and enhance the continuity and accuracy of the inversion results.

[0070] 3. By combining minimum-maximum concavity and higher-order total variation constraints, wave impedance can be reconstructed more accurately in structurally complex regions, significantly improving the resolution and accuracy of the inversion results.

[0071] 4. Since the minimum-maximum concave high-order total variation constraint can effectively suppress the step effect and reduce artifacts in the boundary region, the inversion results are more realistic and reliable.

[0072] 5. This invention significantly reduces the error in the inversion results, thereby improving the accuracy of the seismic impedance inversion results.

[0073] 6. In regions with complex impedance structures, this invention can significantly improve resolution, making the inversion results clearer and more accurate under complex geological conditions. Attached Figure Description

[0074] Figure 1 This is a flowchart of the seismic wave impedance inversion method according to an embodiment of the present invention;

[0075] Figure 2 This is a flowchart illustrating the signal processing of an embodiment of the present invention.

[0076] Figure 3 The following are the impedance inversion results and error profiles of the Marmousi2 model using the traditional total variation algorithm in this embodiment of the invention, where (a) is the impedance profile and (b) is the error profile.

[0077] Figure 4 The images show the impedance inversion results and error profiles of the Marmousi2 model proposed in this embodiment of the invention, where (a) is the impedance profile and (b) is the error profile. Detailed Implementation

[0078] To make the objectives, technical solutions, and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the accompanying drawings and examples.

[0079] like Figure 1 As shown,

[0080] A seismic impedance inversion method using minimum-maximum concavity higher-order total variational constraints includes the following steps:

[0081] S1. Preprocess the seismic data to obtain the initial data for seismic inversion;

[0082] S101. Stack the pre-stack seismic data according to the incident angle to obtain seismic data with limited angles or post-stack seismic data (hereinafter referred to as seismic data). The seismic data is denoted as S.

[0083] S102. Obtain the impedance data corresponding to the target formation section using well logging data. The impedance at a single well logging point is denoted as Z0.

[0084] S103. Using the seismic data from the wellbore and the target layer impedance at the well logging point, obtain the target layer seismic wavelet data w, w = [w1, w2, ..., w p ], where the length of the seismic wavelet data is p, w i (i = 1, 2, ..., p) represents the wavelet data value at the i-th sampling time;

[0085] S104. Based on the target formation stratigraphic data, extrapolate and interpolate the logging data according to the stratigraphic position and perform low-pass filtering to obtain the initial model profile of the target stratigraphic impedance.

[0086] S2. Construct the seismic impedance inversion objective function constrained by minimum and maximum concavity high-order total variational constraints;

[0087] S201. Calculate the initial natural logarithm of the wave impedance L0. The formula for calculating the natural logarithm is:

[0088] L0 = ln(Z0)

[0089] Where Z0 represents the initial wave impedance, L0 represents the natural logarithmic wave impedance, and ln represents the natural logarithmic operation formula;

[0090] S202. Construct a linear calculation model for the reflection coefficient R and its natural logarithm of wave impedance L, the calculation formula is as follows:

[0091] R = DL

[0092] In the above formula, R is the reflection coefficient, L is the natural logarithm of the wave impedance, and D is the difference matrix.

[0093] L = ln(Z)

[0094]

[0095] m represents the total number of sampling points for the logarithm L of the wave impedance;

[0096] S203. Construct the minimum and maximum concavity higher-order total variation constraints. The minimum and maximum concavity higher-order total variation constraints are calculated as follows:

[0097]

[0098] Among them, D xx D represents the second-order difference in the transverse direction. yy Let D represent the longitudinal second-order difference, γ represent the minimum-maximum concavity smoothing constraint parameter (whose value needs to be estimated in advance), ν represent the Moro envelope function factor, ||·||2 represents the L2 norm calculation formula, and ||·||1 represents the L1 norm calculation formula, where D xx L and LD yy The calculation formulas are as follows:

[0099] (D xx L) s,t =2L s,t -L s-1,t -L s+1,t

[0100] (LD yy ) s,t =2L s,t -L s,t-1 -L s,t+1

[0101] S204. Construct a seismic impedance forward model using minimum-maximum concavity higher-order total variational constraints:

[0102]

[0103] In the above formula, S represents seismic data, W represents the seismic wavelet convolution matrix, which is composed of the seismic wavelet data w described in step S103. The calculation method of the minimum-maximum concavity higher-order total variation constraint term is the same as in step S203. λ represents the sparse constraint weight coefficient, the value of which needs to be specified in advance. α represents the initial model constraint weight coefficient, the value of which needs to be specified in advance. The structure of the wavelet convolution matrix W is as follows:

[0104]

[0105] p represents the length of the seismic wavelet w.

[0106] S205. Construct an augmented Lagrangian function in the seismic wave impedance forward model to transform the constrained optimization problem into an unconstrained optimization problem:

[0107]

[0108] In the above formula, R xx and R yy Represents the Lagrange multiplier, C xx and C yy For the dual term, β is a parameter greater than 0.

[0109] S3. Iteratively update the logarithm of wave impedance using the alternating direction multiplier method;

[0110] S301. Use the logarithm of the initial model wave impedance L0 as the natural logarithm of the wave impedance to iteratively calculate the initial value L. 0 Take the initial value R of the Lagrange multipliers. xx =R yy =0 and the initial value of the dual term C xx =C yy =0;

[0111] S302. Update the natural logarithm L of the wave impedance using the solution to the Sylvester equation. k+1 The updated calculation formula is as follows:

[0112]

[0113] Among them, L k+1 Let represent the (k+1)th iteration of the logarithm of wave impedance, vec denotes the transformation of the matrix into a column vector, E represents the identity matrix, and T represents the transpose of the matrix. This represents the Kronecker product operation, where A, B, and C are represented as follows:

[0114]

[0115] S303. Update the Lagrange multiplier R xx and R yy The updated calculation formula is as follows:

[0116]

[0117] in, and This represents the result of the Lagrange multipliers in the (k+1)th iteration, where sign represents the sign function. and The calculation formula is:

[0118]

[0119] Among them, L k+1 This represents the (k+1)th iteration of the logarithm of the wave impedance. and This represents the result of the Lagrange multipliers in the k-th iteration. and This represents the result of the dual term in the k-th iteration;

[0120] S304. Update dual term C xx and C yy The updated calculation formula is as follows:

[0121]

[0122] in and This represents the result of the dual term in the (k+1)th iteration. and Let L represent the result of the dual term in the k-th iteration. k+1 This represents the (k+1)th iteration of the logarithm of the wave impedance. and This represents the result of the Lagrange multiplier subterms in the (k+1)th iteration.

[0123] S4. Calculate the inversion results of seismic wave impedance based on the conditions.

[0124] S401. Update the logarithm L of the seismic impedance in the (k+1)th iteration. k+1 ;

[0125] S402. If ||L k+1 -L k ||2 / ||L k ||2>ε, where ε=5×10 -5 If the condition is met, return to step S302 for a loop; otherwise, calculate the wave impedance Z. The formula for calculating wave impedance Z is:

[0126] Z = exp(L k+1 )

[0127] Here, exp represents the natural exponent operation.

[0128] like Figure 2 As shown, a seismic impedance inversion signal processing flow using minimum-maximum concavity higher-order total variational constraints includes the following steps:

[0129] Obtain the initial value of the wave impedance and calculate its natural logarithm;

[0130] Construct a wave impedance forward model with minimum and maximum concavity higher-order total variation constraints, including data fidelity terms, minimum and maximum concavity higher-order total variation constraints, and initial model constraints;

[0131] The logarithm of the initial model wave impedance is used as the initial value for iterative calculation of the natural logarithm of the wave impedance, and the initial values ​​of the Lagrange multipliers and dual terms are set to 0.

[0132] Update the wave impedance logarithm, Lagrange multipliers, and dual terms using the alternating direction multiplier method;

[0133] Determine whether the updated logarithm of the wave impedance satisfies the loop condition. If so, use the updated logarithm of the wave impedance as the initial value; otherwise, calculate the wave impedance based on the logarithm of the wave impedance.

[0134] To verify the effectiveness of the method proposed in this invention, it was applied to the impedance inversion calculation of the Marmousi2 model and compared with the inversion results of the traditional Total Variation (ATV) method. The inversion results are as follows: Figure 3 and Figure 4 As shown, where Figure 3 (a) and Figure 3 (b) shows the impedance inversion results and impedance error profile of the Marmousi2 model based on the traditional Total Variation (ATV) method. Figure 4 (a) and Figure 4 (b) shows the impedance inversion results and impedance error profile of the proposed minimum maximum concavity higher-order total variation constraint algorithm for the marmousi2 model. As can be seen from the figure, the inversion result error of the proposed method is significantly lower than that of the traditional method, the resolution is significantly improved in regions with complex impedance structures, and the boundary artifacts are suppressed.

[0135] In another embodiment of the present invention, a seismic impedance inversion system is provided, which can be used to implement the above-described seismic impedance inversion method, specifically including:

[0136] Data preprocessing module: Used to preprocess seismic data, including acquiring post-stack seismic record data, well logging data, seismic wavelets, and initial seismic inversion models.

[0137] Inversion objective function construction module: Creates a seismic wave impedance inversion calculation objective function based on data fidelity terms, minimum and maximum concavity high-order total variation constraints, and initial model constraints.

[0138] Inversion Iterative Update Module: The natural logarithm of wave impedance is updated based on the alternating direction multiplier method and the objective function of seismic wave impedance inversion, resulting in the updated natural logarithm of wave impedance.

[0139] Inversion Result Judgment Module: Determines whether the values ​​of the natural logarithm of wave impedance before and after the update meet the iterative update conditions. If so, the updated natural logarithm of wave impedance, Lagrange multipliers, and dual terms are used as initial values, and the module is transferred to the inversion iterative update module. If not, the updated natural logarithm of wave impedance is subjected to exponential operation to obtain the inversion wave impedance result.

[0140] Results Output Module: Outputs and stores the final inverted wave impedance results, including generating visualizations and detailed reports, ensuring that all output data and charts meet the expected format and accuracy requirements.

[0141] In another embodiment of the present invention, a terminal device is provided, comprising a processor and a memory. The memory stores a computer program, which includes program instructions. The processor executes the program instructions stored in the computer storage medium. The processor may be a Central Processing Unit (CPU), or other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. It is the computing and control core of the terminal, suitable for implementing one or more instructions, specifically suitable for loading and executing one or more instructions to achieve a corresponding method flow or corresponding function. The processor described in this embodiment of the present invention can be used for the operation of a seismic impedance inversion method.

[0142] In another embodiment of the present invention, a storage medium is provided, specifically a computer-readable storage medium (Memory). This computer-readable storage medium is a memory device in a terminal device used to store programs and data. It is understood that the computer-readable storage medium here can include both the built-in storage medium in the terminal device and extended storage media supported by the terminal device. The computer-readable storage medium provides storage space that stores the terminal's operating system. Furthermore, this storage space also stores one or more instructions suitable for loading and execution by a processor. These instructions can be one or more computer programs (including program code). It should be noted that the computer-readable storage medium here can be high-speed RAM or non-volatile memory, such as at least one disk storage device.

[0143] One or more instructions stored in a computer-readable storage medium can be loaded and executed by a processor to implement the corresponding steps of the seismic wave impedance inversion method in the above embodiments; one or more instructions in the computer-readable storage medium are loaded and executed by a processor.

[0144] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0145] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0146] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0147] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0148] Those skilled in the art will recognize that the embodiments described herein are intended to help the reader understand the implementation methods of the present invention, and should be understood that the scope of protection of the present invention is not limited to such specific statements and embodiments. Those skilled in the art can make various other specific modifications and combinations based on the technical teachings disclosed in this invention without departing from the spirit of the invention, and these modifications and combinations are still within the scope of protection of the present invention.

Claims

1. A seismic wave impedance inversion method with minimum-maximum concavity higher-order total variational constraints, characterized in that, Includes the following steps: S1. Preprocess the seismic data to obtain the initial data for seismic inversion, including: finite angle or post-stack seismic data, impedance data corresponding to the target stratigraphic segment, seismic wavelet data of the target layer, and the initial model for seismic inversion; S2. Construct the seismic impedance inversion objective function with minimum and maximum concavity higher-order total variation constraints, including: constructing a linear relationship model between reflection coefficient and natural logarithm of wave impedance based on data fidelity terms; constructing sparse constraint terms based on minimum and maximum concavity higher-order total variation constraints; and constructing the seismic impedance inversion objective function based on initial model constraint terms. S3. Iteratively update the logarithm of wave impedance using the alternating direction multiplier method, calculate and update the natural logarithm of wave impedance, and update the Lagrange multiplier terms and dual terms; S4. Calculate the inversion result of seismic wave impedance according to the conditions, and determine whether the values ​​before and after updating the natural logarithm of wave impedance meet the iterative update conditions. If so, take the updated natural logarithm of wave impedance, Lagrange multiplier term and dual term as the initial values ​​and transfer to the inversion iterative update module. If not, perform exponential operation on the updated natural logarithm of wave impedance to obtain the inversion wave impedance result. S5. Output and store the final inverted wave impedance results for easy subsequent analysis and application.

2. The seismic wave impedance inversion method according to claim 1, characterized in that: Step S1 includes the following sub-steps: S101. Stack pre-stack seismic data according to the incident angle to obtain seismic data with limited angles or post-stack seismic data, referred to as seismic data, and the seismic data is denoted as S; S102. Obtain the impedance data corresponding to the target formation section using well logging data. The impedance at a single well logging point is denoted as Z0. S103. Using the seismic data from the wellbore and the target layer impedance at the well logging point, obtain the target layer seismic wavelet data w, w = [w1, w2, ..., w p ], where the length of the seismic wavelet data is p, w i (i = 1, 2, ..., p) represents the wavelet data value at the i-th sampling time; S104. Based on the target formation stratigraphic data, perform low-pass filtering on the well logging data and extrapolate the values ​​according to the stratigraphic position to obtain the initial model profile of the wave impedance of the target stratigraphic segment.

3. The seismic wave impedance inversion method according to claim 1, characterized in that: Step S2 includes the following sub-steps: S201. Calculate the initial natural logarithm of the wave impedance L0. The formula for calculating the natural logarithm is: L0 = ln(Z0) Where Z0 represents the initial wave impedance, L0 represents the natural logarithmic wave impedance, and ln represents the natural logarithmic operation formula; S202. Construct a linear calculation model for the reflection coefficient R and its natural logarithm of wave impedance L, the calculation formula is as follows: R = DL In the above formula, R is the reflection coefficient, L is the natural logarithm of the wave impedance, and D is the difference matrix. L = ln(Z) m represents the total number of sampling points for the logarithm L of the wave impedance; S203. Construct the minimum and maximum concavity higher-order total variation constraints. The minimum and maximum concavity higher-order total variation constraints are calculated as follows: Among them, D xx D represents the second-order difference in the transverse direction. yy Let D represent the longitudinal second-order difference, γ represent the minimum-maximum concavity smoothing constraint parameter (whose value needs to be estimated in advance), ν represent the Moro envelope function factor, ||·||2 represents the L2 norm calculation formula, and ||·||1 represents the L1 norm calculation formula, where D xx L and LD yy The calculation formulas are as follows: (D xx L) s,t <2L s,t -L s-1,t -L s+1,t (LD yy ) s,t <2L s,t -L s,t-1 -L s,t+1 S204. Construct a seismic impedance forward model using minimum-maximum concavity higher-order total variational constraints: In the above formula, S represents seismic data, W represents the seismic wavelet convolution matrix, which is composed of the seismic wavelet data w described in step S103. The calculation method of the minimum-maximum concavity higher-order total variation constraint term is the same as in step S203. λ represents the sparse constraint weight coefficient, the value of which needs to be specified in advance. α represents the initial model constraint weight coefficient, the value of which needs to be specified in advance. The structure of the wavelet convolution matrix W is as follows: p represents the length of the seismic wavelet w; S205. Construct an augmented Lagrangian function in the seismic wave impedance forward model to transform the constrained optimization problem into an unconstrained optimization problem: In the above formula, R xx and R yy Represents the Lagrange multiplier, C xx and C yy For the dual term, β is a parameter greater than 0.

4. The seismic wave impedance inversion method according to claim 1, characterized in that: Step S3 includes the following sub-steps: S301. Use the logarithm of the initial model wave impedance L0 as the natural logarithm of the wave impedance to iteratively calculate the initial value L. 0 Take the initial value R of the Lagrange multipliers. xx =R yy =0 and the initial value of the dual term C xx =C yy =0; S302. Update the natural logarithm L of the wave impedance using the solution to the Sylvester equation. k+1 The updated calculation formula is as follows: Among them, L k+1 Let represent the (k+1)th iteration of the logarithm of wave impedance, vec denotes the transformation of the matrix into a column vector, E represents the identity matrix, and T represents the transpose of the matrix. This represents the Kronecker product operation, where A, B, and C are represented as follows: S303. Update the Lagrange multiplier R xx and R yy The updated calculation formula is as follows: in, and This represents the result of the Lagrange multipliers in the (k+1)th iteration, where sign represents the sign function. and The calculation formula is: Among them, L k+1 This represents the (k+1)th iteration of the logarithm of the wave impedance. and This represents the result of the Lagrange multipliers in the k-th iteration. and This represents the result of the dual term in the k-th iteration; S304. Update dual term C xx and C yy The updated calculation formula is as follows: in and This represents the result of the dual term in the (k+1)th iteration. and Let L represent the result of the dual term in the k-th iteration. k+1 This represents the (k+1)th iteration of the logarithm of the wave impedance. and This represents the result of the Lagrange multiplier subterms in the (k+1)th iteration.

5. The seismic wave impedance inversion method according to claim 4, characterized in that: Step S4 includes the following sub-steps: S401. Update the logarithm L of the seismic impedance in the (k+1)th iteration. k+1 ; S402. If ||L k+1 -L k ||2 / ||L k ||2>ε, where ε=5×10 -5 If the condition is met, return to step S302 for a loop; otherwise, calculate the wave impedance Z. The formula for calculating wave impedance Z is: Z=exp(L k+1 ) Here, exp represents the natural exponent operation.

6. A seismic wave impedance inversion system, characterized in that: This system can be used to implement the seismic wave impedance inversion method according to any one of claims 1 to 5, specifically including: Data preprocessing module: used to preprocess seismic data, including acquiring post-stack seismic record data, well logging data, seismic wavelets, and initial seismic inversion models; Inversion objective function construction module: Based on the data fidelity term, the minimum and maximum concavity high-order total variation constraint term, and the initial model constraint term, a seismic wave impedance inversion calculation objective function is created; Inversion Iterative Update Module: The natural logarithm of wave impedance is updated based on the alternating direction multiplier method and the objective function of seismic wave impedance inversion, resulting in the updated natural logarithm of wave impedance; Inversion Result Judgment Module: Determine whether the values ​​of the natural logarithm of wave impedance before and after the update meet the iterative update conditions. If so, use the updated natural logarithm of wave impedance, Lagrange multipliers and dual terms as initial values ​​and transfer to the inversion iterative update module. If not, perform exponential operation on the updated natural logarithm of wave impedance to obtain the inversion wave impedance result. Results Output Module: Outputs and stores the final inverted wave impedance results, including generating visualizations and detailed reports, ensuring that all output data and charts meet the expected format and accuracy requirements.

7. A computer device, characterized in that: It includes a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor, when executing the program, implements the seismic impedance inversion method according to any one of claims 1 to 5.

8. A computer-readable storage medium, characterized in that: It stores a computer program that, when executed by a processor, implements the seismic impedance inversion method according to any one of claims 1 to 5.

Images