Inversion method and system based on hybrid regularization model, electronic equipment and storage medium

Abstract

The invention provides an inversion method and system based on a hybrid regularization model, electronic equipment and a storage medium. The method comprises the following steps: S1, inputting seismic data, and constructing a target function; s2, judging whether the number of iterations of the target function reaches the maximum value Kmax and whether the error of two times of iterations reaches a small value or not; s3, calculating a pseudo gradient # imgabs0 # of the target function according to a judgment result of the step 2; S4, calculating a descending direction of the pseudo gradient # imgabs1 # by utilizing L-BFGS loop iteration, and mapping the same quadrant for correction; and S5, according to the descending direction of the pseudo gradient # imgabs2, calculating a new inversion result through a linear search strategy. According to the method, through mixed regularization constraint inversion, the development condition of the stratum is better described, and the OWL-QN algorithm is introduced to accelerate the convergence of the algorithm and improve the inversion precision.

Description

Technical Field

[0001] The present invention belongs to the technical field of seismic inversion, and particularly relates to an inversion method, system, electronic device and storage medium based on a hybrid regularization model. Background Technique

[0002] Currently, seismic inversion methods can effectively identify the lithology and hydrocarbon-bearing characteristics of reservoirs. This technology is widely used to estimate formation elastic parameters. However, inverse problems are often accompanied by ill-posedness, mainly manifested as the problems of non-uniqueness and sensitivity to noise. To address these problems, regularization techniques are often used to obtain stable approximate solutions.

[0003] In the mid-20th century, mathematician Tikhonov conducted research on inverse problems and proposed the Tikhonov regularization method, which is a method that uses a norm to describe the regularization term of the inversion to constrain the inversion. Therefore, the Tikhonov regularization method cannot obtain boundary features. The TV regularization method can adapt to and describe such boundary features. As a boundary-preserving regularization method, by using the L1 norm constraint of the gradient of the model parameters to be inverted, it is expected to effectively identify the geological environment with abrupt changes in formation parameter attributes. Then, in the context of complex geological conditions, in order to obtain and protect formation boundary features and retain smooth information, a hybrid constraint strategy is often adopted to better characterize formation information, and the expression can be written as:

[0004]

[0005] where G is the forward matrix, m is the vector, d is the target value, and L is the difference matrix. This is a convex optimization problem, where α, β > 0 are regularization parameters used to adjust the quadratic norm term and the regularization constraint term of the data. To find the solution in Equation (1), the gradient method or the Newton method is usually used to solve it, so that the objective function can better converge to the optimal solution.

[0006] The Newton method estimates the minimum point of a twice-differentiable real function. Its outstanding advantage is fast convergence. However, applying the Newton method requires calculating second-order partial derivatives, and the Hessian matrix of the objective function may be non-positive definite. To overcome this defect, the quasi-Newton method was proposed. Its basic idea is to approximate the inverse matrix of the Hessian matrix in the Newton method with a matrix that does not contain second-order derivatives. For its specific form, in 1970, four mathematicians Broyden, Fletcher, Goldfard, and Shanno named the BFGS algorithm with their names. Then, for the case where the scale of the optimization problem to be solved is very large, the limited-memory BFGS algorithm, that is, the L-BFGS algorithm, is proposed for optimization calculation.

[0007] In the above method, the inverse problem is often accompanied by ill-posedness, mainly manifested as the problem of multiple solutions and sensitivity to noise.

[0008] The above technical problems need to be solved urgently. Summary of the Invention

[0009] To solve the above technical problems, the present invention proposes a technical solution of an inversion method based on a hybrid regularization model to improve the accuracy and efficiency of reservoir inversion and solve the above technical problems.

[0010] The first aspect of the present invention discloses an inversion method based on a hybrid regularization model, and the method includes:

[0011] Step S1: Input seismic data and construct an objective function;

[0012] Step S2: Determine whether the number of iterations of the objective function reaches the maximum value Kmax and whether the error between the previous and the current iteration reaches a small value;

[0013] Step S3: According to the judgment result of step 2, calculate the pseudo-gradient ◇g(m i );

[0014] Step S4: Use the L-BFGS cyclic iteration to calculate the descent direction of the pseudo-gradient ◇g(m i ) and map it to the same quadrant for correction;

[0015] Step S5: According to the descent direction of the pseudo-gradient ◇g(m i ), calculate a new inversion result through a linear search strategy.

[0016] According to the method of the first aspect of the present invention, in step S1, when the seismic wave is vertically incident, the reflection coefficient of adjacent formation interfaces can be represented by the acoustic impedance of the formation, and the formation acoustic impedance is continuous. Therefore, the reflection coefficient of the formation is expressed as:

[0017]

[0018] For a single seismic trace, the specific formula is:

[0019] d = Gm + n

[0020] where d ∈ R N represents single-trace seismic data, which is a vector about the time series t; G is the forward operator, composed of W*L, W = R N×N is the convolution operator composed of the source wavelet; L is the finite difference operator, m ∈ R N represents the elastic parameter vector of a single-trace formation; n ∈ R N represents the random noise vector;

[0021] Therefore, the objective function is as follows:

[0022]

[0023] According to the method of the first aspect of the present invention, in step S2, when the number of iterations of the objective function reaches the maximum value Kmax and the error between two consecutive times reaches a small value, the iteration is stopped; otherwise, the iteration continues.

[0024] According to the method of the first aspect of the present invention, in step S3, when the iteration in step S2 stops, the pseudo-gradient ◇g(m i ) of the objective function is calculated, and the calculation method is as follows:

[0025]

[0026] where g i (m) is the gradient of the differentiable L2-norm term in the objective function, sign(*) represents the vector sign function, ξ i represents the i-th element after the L1-norm term is vectorized, and β is the regularization parameter of the total variation term.

[0027] According to the method of the first aspect of the present invention, in step S4, for solving large-scale non-linear optimization, the inverse of the approximate Hessian matrix is calculated by L-BFGS, so as to obtain the descent direction of the pseudo-gradient, that is:

[0028] d k =-H k ◇g(m k )

[0029] where H k is the inverse of the Hessian matrix of the objective function calculated by L-BFGS; ◇g(m i ) is the pseudo-gradient of the objective function; k is the number of iterations;

[0030] According to the descent direction of the pseudo-gradient, the vector before and after can be updated to the same quadrant, and the orthogonal function is defined as:

[0031]

[0032] where sign(*) represents the vector sign function; x i represents the formation parameter result of the first vector in the i-th layer;

[0033] To ensure that the two variables fall into the same quadrant, the descent direction of the new pseudo-gradient is:

[0034] p k =π(d k, -◇g(m k ))。

[0035] According to the method of the first aspect of the present invention, in the step S5, the step size α of the solution result k adopts a linear search strategy as follows:

[0036]

[0037] where when when the step size α k is the step size calculated by adopting the linear search strategy; represents the formation parameter of the i-th layer at the k-th iteration, represents one of the vectors of the orthogonal function, and sign(·) represents the sign of the vector.

[0038] According to the method of the first aspect of the present invention, in the step S5, when the number of iterations k reaches the maximum or the difference between the results of two iterations is within a certain range, the algorithm is terminated, and the corresponding inversion result is output.

[0039] The second aspect of the present invention discloses an inversion system based on a hybrid regularization model, and the system includes:

[0040] A first processing module, configured to input seismic data and construct an objective function;

[0041] A second processing module, configured to judge whether the number of iterations of the objective function reaches the maximum value Kmax and whether the error between the previous and the next time reaches a small value;

[0042] A third processing module, configured to calculate the pseudo-gradient ◇g(m i ) of the objective function according to the judgment result of step 2;

[0043] A fourth processing module, configured to use the L-BFGS cyclic iteration to calculate the descent direction of the pseudo-gradient ◇g(m i ) and map and correct it in the same quadrant;

[0044] A fifth processing module, configured to calculate a new inversion result through a linear search strategy according to the descent direction of the pseudo-gradient ◇g(m i ).

[0045] According to the system of the second aspect of the present invention, the first processing module is specifically configured that when the seismic wave is vertically incident, the reflection coefficient of the adjacent formation interface can be represented by the acoustic impedance of the formation, and the formation acoustic impedance is continuous. Therefore, the reflection coefficient of the formation is expressed as:

[0046]

[0047] For a single seismic trace, the specific formula is:

[0048] d = Gm + n

[0049] where d ∈ R N represents the single-trace seismic data, which is a vector with respect to the time series t; G is the forward operator, composed of W * L, and W = R N×N is the convolution operator composed of the source wavelet; L is the finite difference operator, and m ∈ R N represents the elastic parameter vector of the single-trace formation; n ∈ R N represents the random noise vector;

[0050] Therefore, the objective function is:

[0051]

[0052] According to the system of the second aspect of the present invention, the second processing module is specifically configured to stop iteration when the number of iterations of the objective function reaches the maximum value Kmax and the error between the previous and the next time reaches a small value; otherwise, continue the iteration.

[0053] According to the system of the second aspect of the present invention, the third processing module is specifically configured to calculate the pseudo-gradient ◇g(m i ) when the iteration stops, and the calculation method is:

[0054]

[0055] where g i (m) is the gradient of the differentiable quadratic term in the objective function, sign(*) represents the vector sign function, ξ i represents the i-th element after vectorizing the L1 norm term, and β is the regularization parameter of the total variation term.

[0056] According to the system of the second aspect of the present invention, the fourth processing module is specifically configured to, for solving large-scale non-linear optimization, calculate the inverse of the approximate Hessian matrix through L-BFGS, so as to obtain the descent direction of the pseudo-gradient, that is:

[0057] d k = -H k ◇g(m k )

[0058] where H k is the inverse of the Hessian matrix of the objective function calculated by L-BFGS; ◇g(m i) is the pseudo-gradient of the objective function; k is the number of iterations;

[0059] According to the descending direction of the pseudo-gradient, if the vector mappings before and after the update are mapped to the same quadrant, the orthogonal function is defined as:

[0060]

[0061] where sign(*) represents the vector sign function; x i represents the formation parameter result of the first vector in the i-th layer;

[0062] To ensure that the two variables fall into the same quadrant, the descending direction of the new pseudo-gradient is:

[0063] p k = π(d k , -◇g(m k ))

[0064] According to the system of the second aspect of the present invention, the fifth processing module is specifically configured to solve the step size α of the result k The strategy adopted for linear search is:

[0065]

[0066] where when when The step size α k is the step size calculated by adopting the strategy of linear search; represents the formation parameter of the i-th layer at the k-th iteration, represents one of the vectors of the orthogonal function, and sign(·) represents the sign of the vector.

[0067] According to the system of the second aspect of the present invention, the fifth processing module is specifically configured to end the algorithm and output the corresponding inversion result when the number of iterations k reaches the maximum or the difference between the results of two iterations is within a certain range.

[0068] The third aspect of the present invention discloses an electronic device. The electronic device includes a memory and a processor. When the processor executes the computer program stored in the memory, the steps in any one of the inversion methods of a basic hybrid regularization model in the first aspect of the present disclosure are implemented.

[0069] The fourth aspect of the present invention discloses a computer-readable storage medium. A computer program is stored on the computer-readable storage medium. When the computer program is executed by the processor, the steps in any one of the inversion methods of a basic hybrid regularization model in the first aspect of the present disclosure are implemented.

[0070] In summary, the solution proposed by the present invention can better characterize the development of the formation through hybrid regularization constrained inversion, and introduce the OWL-QN algorithm to accelerate the convergence of the algorithm and improve the accuracy of inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[0071] In order to more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the following will briefly introduce the drawings required for the description of the specific embodiments or the prior art. Obviously, the drawings in the following description are some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can also be obtained based on these drawings.

[0072] Figure 1 FIG. 9 is a flowchart of an inversion method based on a hybrid regularization model according to an embodiment of the present invention;

[0073] Figure 2 FIG. 13 is a flowchart of a specific embodiment of an inversion method based on a hybrid regularization model according to an embodiment of the present invention;

[0074] FIG. 3(a) shows the inversion result of the block model under the condition of no noise;

[0075] FIG. 3(b) shows the inversion result of the block model under the condition that the signal-to-noise ratio is 2;

[0076] FIG. 4(a) shows the inversion result of the curve model under the condition of no noise;

[0077] FIG. 4(b) shows the inversion result of the curve model under the condition that the signal-to-noise ratio is 2;

[0078] FIG. 5(a) is the initial model of the Marmousi model;

[0079] FIG. 5(b) is the true model of the Marmousi model;

[0080] Figure 6 FIG. 35 shows the inversion result of the impedance attribute under the hybrid regularization constraint according to an embodiment of the present invention;

[0081] Figure 7 FIG. 39 shows the convergence comparison and error analysis of the hybrid regularization and the total variation regularization according to an embodiment of the present invention;

[0082] Figure 8 FIG. 43 is a structural diagram of an inversion system based on a hybrid regularization model according to an embodiment of the present invention;

[0083] Figure 9 FIG. 47 is a structural diagram of an electronic device according to an embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0084] To make the objectives, technical solutions and advantages of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Apparently, the described embodiments are some, but not all, of the embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0085] Embodiment 1:

[0086] The first aspect of the present invention discloses an inversion method based on a hybrid regularization model. Figure 1 As shown in the flowchart of an inversion method based on a hybrid regularization model according to an embodiment of the present invention, Figure 1 as shown, the method includes:

[0087] Step S1: Input seismic data and construct an objective function;

[0088] Step S2: Determine whether the number of iterations of the objective function reaches the maximum value Kmax and whether the error between two consecutive times reaches a small value;

[0089] Step S3: According to the judgment result of Step 2, calculate the pseudo-gradient ◇g(m i );

[0090] Step S4: Use the L-BFGS cyclic iteration to calculate the descent direction of the pseudo-gradient ◇g(m i ) and map it to the same quadrant for correction;

[0091] Step S5: According to the descent direction of the pseudo-gradient ◇g(m i ), calculate a new inversion result through a linear search strategy.

[0092] In some embodiments, in Step S1, when the seismic wave is vertically incident, the reflection coefficient of adjacent formation interfaces can be represented by the acoustic impedance of the formation. The formation acoustic impedance is continuous. Therefore, the reflection coefficient of the formation is expressed as:

[0093]

[0094] For a single seismic trace, the specific formula is:

[0095] d = Gm + n

[0096] where d ∈ R N represents single-trace seismic data and is a vector with respect to the time series t; G is the forward operator, composed of W*L, and W = R N×N is the convolution operator composed of the source wavelet; L is the finite difference operator, and m ∈ RN denotes the elastic parameter vector of a single-layer formation; n ∈ R N denotes the random noise vector;

[0097] Therefore, the objective function is:

[0098]

[0099] In some embodiments, in step S2, when the number of iterations of the objective function reaches the maximum value Kmax and the error between two consecutive times reaches a small value, the iteration stops; otherwise, the iteration continues.

[0100] In some embodiments, in step S3, when the iteration in step S2 stops, calculate the pseudo-gradient ◇g(m i ), and the calculation method is:

[0101]

[0102] where g i (m) is the gradient of the differentiable quadratic term in the objective function, sign(*) represents the vector sign function, ξ i represents the i-th element after vectorizing the L1 norm term, and β is the regularization parameter of the total variation term.

[0103] In some embodiments, in step S4, for solving large-scale non-linear optimization, calculate the inverse of the approximate Hessian matrix through L-BFGS, so as to obtain the descent direction of the pseudo-gradient, that is:

[0104] d k =-H k ◇g(m k )

[0105] where H k is the inverse of the Hessian matrix of the objective function calculated by L-BFGS; ◇g(m i ) is the pseudo-gradient of the objective function; k is the number of iterations;

[0106] According to the descent direction of the pseudo-gradient, the vectors before and after can be updated to the same quadrant, and the orthogonal function is defined as:

[0107]

[0108] where sign(*) represents the vector sign function; x i represents the formation parameter result of the first vector in the i-th layer; to ensure that the two variables fall into the same quadrant, the descent direction of the new pseudo-gradient is:

[0109] pk = π(d k , -◇g(m k ))。

[0110] In some embodiments, in the step S5, the step size α of the solution result k adopts a linear search strategy as follows:

[0111]

[0112] where when when the step size α k is the step size calculated by adopting the linear search strategy; represents the formation parameter of the i-th layer at the k-th iteration, represents one of the vectors of the orthogonal function, and sign(·) represents the sign of the vector.

[0113] In some embodiments, in the step S5, when the number of iterations k reaches the maximum or the difference between the results of two iterations is within a certain range, the algorithm ends and the corresponding inversion result is output.

[0114] In summary, the present invention can improve the accuracy, anti-noise performance and calculation efficiency of inversion, and uses a hybrid regularization term to constrain inversion. Due to the characteristics of model constraints, the present invention can better protect the formation boundary and internal characteristics of the formation, and effectively depict the development of the formation.

[0115] Embodiment 2:

[0116] For the case where the scale of the optimization problem to be solved is very large, the prior memory BFGS algorithm, that is, the L-BFGS algorithm, is proposed in the prior art for optimization calculation. On this basis, the orthogonal limited memory quasi-Newton algorithm (OWL-QN) algorithm is developed, that is, the L-BFGS algorithm is updated under the assumption of the variable quadrant, so that the front and rear variables are in the same quadrant. The present invention uses this algorithm for solution; as Figure 2 shown, a solution process for elastic parameter inversion under hybrid regularization constraints is provided, which can consider the seismic inversion method under hybrid regularization constraints, better retain a boundary feature, as well as the internal smoothing result, and highlight a change feature in the vertical direction, thereby improving the inversion accuracy.

[0117] When the seismic wave is vertically incident, the reflection coefficient of the adjacent formation interface can be represented by the acoustic impedance of the formation. Foster believes that the formation acoustic impedance is continuous, so the reflection coefficient of the formation can be expressed as:

[0118]

[0119] For a single seismic trace, it can be written as:

[0120] d = Gm + n (2)

[0121] where d ∈ R N represents the single-trace seismic data, which is a vector with respect to the time series t; G is the forward operator, composed of W * L, where W = R N×N is the convolution operator composed of the source wavelet; L is the finite-difference operator, m ∈ R N represents the elastic parameter vector of the single-trace formation; n ∈ R N represents the random noise vector.

[0122] Then the objective function can be constructed, that is:

[0123]

[0124] For a large-scale non-linear optimization problem like (4), the L-BFGS method can be used to approximately calculate the inverse of the Hessian matrix, which can not only improve the inversion speed but also avoid the memory consumption caused by storing large-scale matrices. For the non-differentiable L1-norm term in the objective functional (1), the OWL-QN solution algorithm can be introduced. By fixing the quadrant, the objective function can be made differentiable. It has only a few differences compared with the standard L-BFGS algorithm, that is, using the pseudo-gradient to replace the gradient, the search direction needs to be consistent with the pseudo-gradient direction to keep the quadrant unchanged, and the next selected point has the same quadrant as the current point. These conditions can enable the algorithm to better find the target solution. The following is the algorithm introduction:

[0125] First, calculate the pseudo-gradient of the objective function. The calculation method is:

[0126]

[0127] where g i (m) is the gradient of the differentiable quadratic term in the objective function, and sign(*) represents the vector sign function. ξ i represents the i-th element after vectorizing the L1-norm term.

[0128] For solving the large-scale non-linear optimization problem, L-BFGS calculates the inverse of the approximate Hessian matrix to find the descent direction, that is:

[0129] d k = -H k ◇g(m k ) (5)

[0130] where H kis the inverse of the Hessian matrix of the objective function calculated by L - BFGS; ◇g(m i ) is the pseudo - gradient of the objective function; k is the number of iterations.

[0131] With the direction of gradient descent, the algorithm proposes to ensure that the vectors before and after the update are mapped to the same quadrant, and the orthogonal function is defined as:

[0132]

[0133] where sign(*) represents the vector sign function.

[0134] To ensure that the two variables fall into the same quadrant.

[0135] Then there is a new gradient descent direction:

[0136] p k = π(d k , -◇g(m k )) (7)

[0137] And a new solution result, where the step size α of the solution result k adopts a linear search strategy:

[0138]

[0139] where, when when the step size α k is the step size calculated by the linear search strategy.

[0140] When the number of iterations k reaches the maximum or the difference between the results of two iterations is very small, the algorithm can be terminated and the corresponding inversion result can be output.

[0141] As Figure 2 shown, according to the above method, the inversion method of the basis - mixed regularization model is given as follows:

[0142] S101: Read the input seismic data d, wavelet data W, and the initial model of impedance data with horizon constraints generated from well - logging data, that is, construct the objective function;

[0143] S102: Judge whether the number of iterations reaches the maximum value Kmax and whether the error between the previous and the current iterations reaches a minimum value;

[0144] S103: If the result of step S102 is yes, calculate the pseudo - gradient ◇g(m i ) according to equation (5);

[0145] S104: Using L-BFGS loop iteration, calculate the pseudo gradient descent direction in (6), and map it to the same quadrant for correction;

[0146] S105: According to the pseudo gradient ◇g(m i ) is descending, and the new inversion result m is calculated by the linear search strategy. k+1 ;

[0147] In order to verify the stability and adaptability of this method, different model data were used for testing. The test results Figure 3(a) - Figure 3(b) As shown. It can be seen that no matter what kind of model the inversion result is under, the inversion method can obtain a relatively stable and accurate inversion solution. The present invention is applied to block models, curve models and common marmousi models. It can be seen that they all show good inversion results, and the inversion effect of mixed regularization is better than the single model constraint result. From the results in Figure 3, it can be seen that in Figure 3 (a), under the conditions of no noise and the same initial model, the result of the mixed regularization constraint (Invmix) is roughly similar to the result of the total variation regularization constraint (InvTV), and both are basically consistent with the true model (True), while the result under the smoothing constraint (LS) has a strong jitter and brings certain errors; in Figure 3 (b), under the conditions of a signal-to-noise ratio of 2 and the same initial model, the mixed regularization constraint (Invmix) is closer to the true solution (True), while the result under the total variation constraint (InvTV) has a very obvious fast feature, and the smoothing response between layers is weak. In contrast, the smoothing constraint (LS) is closer to the initial model. Then the curve model is inverted, as shown Figure 4(a) - Figure 4(b) As shown, it can be seen that in the case of no noise and the same initial model as Figure 4(a), the result of the mixed regularization constraint (Invmix) is roughly similar to the result of the total variation regularization (InvTV), and both are basically consistent with the true model (True), while the inversion result (LS) under the smoothness constraint will have some deviations in the deep layer. Then for the case of a signal-to-noise ratio of 2 and the same initial model as Figure 4(b), the result of the mixed regularization constraint (Invmix) is closer to the true model (True), while the other two models are relatively deviated. The marmousi model is further tested. Figure 5(a) represents the true model, and Figure 5(b) represents the initial model. In order to verify the lateral continuity and the inversion effect of the profile, it is inverted using the mixed regularization constraint under the condition of a signal-to-noise ratio of 5. The inversion results are shown in the figure below. Figure 6 As shown in the figure, we can see that a good inversion result is obtained, which can fully describe the real model as a whole. Figure 7 We can further see the convergence of the algorithm and the accuracy of the inversion.

[0148] In the present invention, a hybrid regularization objective functional is established, and the orthogonal limited-memory quasi-Newton algorithm (OWL-QN) is adopted for iterative calculation in the solution idea. The introduction of this algorithm not only improves the calculation efficiency of the inversion but also saves memory, making the inversion result closer to the true solution. Different inversion models are used for testing and a reasonable solution is proposed. In the model test, noisy or noise-free data is used for testing, and the elastic parameters of the formation are effectively estimated through the hybrid regularization constraint inversion method, and the comparison with the true formation is carried out. Finally, a trial calculation is performed on the target area to obtain a stable and reasonable inversion result.

[0149] The second aspect of the present invention discloses an inversion system based on a hybrid regularization model. Figure 8 As shown in the structure diagram of an inversion system based on a hybrid regularization model according to an embodiment of the present invention; Figure 8 As shown, the system 100 includes:

[0150] A first processing module 101, configured to input seismic data and construct an objective function;

[0151] A second processing module 102, configured to determine whether the number of iterations of the objective function reaches the maximum value Kmax and whether the error between two consecutive times reaches a small value;

[0152] A third processing module 103, configured to calculate the pseudo-gradient ◇g(m i ) of the objective function according to the judgment result of step 2;

[0153] A fourth processing module 104, configured to use L-BFGS cyclic iteration to calculate the descent direction of the pseudo-gradient ◇g(m i ) and map and correct it in the same quadrant;

[0154] A fifth processing module 105, configured to calculate a new inversion result through a linear search strategy according to the descent direction of the pseudo-gradient ◇g(m i ).

[0155] According to the system of the second aspect of the present invention, the first processing module 101 is specifically configured to, when the seismic wave is vertically incident, the reflection coefficient of adjacent formation interfaces can be represented by the acoustic impedance of the formation, and the formation acoustic impedance is continuous. Therefore, the reflection coefficient of the formation is expressed as:

[0156]

[0157] For a single seismic trace, the specific formula is:

[0158] d = Gm + n

[0159] where d ∈ R Ndenotes single - trace seismic data, which is a vector with respect to the time series t; G is the forward operator, composed of W*L, where W = R N×N is the convolution operator composed of the source wavelet; L is the finite - difference operator, and m ∈ R N denotes the elastic - parameter vector of the single - trace formation; n ∈ R N denotes the random - noise vector;

[0160] Therefore, the objective function is:

[0161]

[0162] According to the system of the second aspect of the present invention, the second processing module 102 is specifically configured to stop iteration when the number of iterations of the objective function reaches the maximum value Kmax and the error between two consecutive times reaches a small value; otherwise, continue the iteration.

[0163] According to the system of the second aspect of the present invention, the third processing module 103 is specifically configured to, when the iteration stops, calculate the pseudo - gradient ◇g(m i ), and the calculation method is:

[0164]

[0165] where g i (m) is the gradient of the differentiable quadratic - norm term in the objective function, sign(*) represents the vector sign function, ξ i represents the i - th element after the L1 - norm term is vectorized, and β is the regularization parameter of the total - variation term.

[0166] According to the system of the second aspect of the present invention, the fourth processing module 104 is specifically configured to, for solving large - scale non - linear optimization, calculate the inverse of the approximate Hessian matrix through L - BFGS, so as to obtain the descent direction of the pseudo - gradient, that is:

[0167] d k =-H k ◇g(m k )

[0168] where H k is the inverse of the Hessian matrix of the objective function calculated by L - BFGS; ◇g(m i ) is the pseudo - gradient of the objective function; k is the number of iterations;

[0169] According to the descent direction of the pseudo - gradient, the vectors before and after can be updated to the same quadrant, and the orthogonal function is defined as:

[0170]

[0171] where, sign(*) represents the vector sign function; x i represents the formation parameter result of the first vector in the i-th layer; to ensure that the two variables fall in the same quadrant, the descent direction of the new pseudo-gradient is:

[0172] p k = π(d k , -◇g(m k ))

[0173] According to the system of the second aspect of the present invention, the fifth processing module 105 is specifically configured to solve the step size α of the result k The strategy adopted for linear search is:

[0174]

[0175] where when when The step size α k is the step size calculated by adopting the strategy of linear search; represents the formation parameter of the i-th layer at the k-th iteration, represents one of the vectors of the orthogonal function, and sign(·) represents the sign of the vector.

[0176] According to the system of the second aspect of the present invention, the fifth processing module 105 is specifically configured to end the algorithm and output the corresponding inversion result when the number of iterations k reaches the maximum or the difference between the results of two iterations is within a certain range.

[0177] The third aspect of the present invention discloses an electronic device. The electronic device includes a memory and a processor. When the processor executes a computer program stored in the memory, the steps in an inversion method of a base hybrid regularization model according to any one of the first aspects disclosed in the present invention are implemented.

[0178] Figure 9 is a structural diagram of an electronic device according to an embodiment of the present invention, as Figure 9As shown in the figure, the electronic device includes a processor, a memory, a communication interface, a display screen, and an input device connected through a system bus. Among them, the processor of the electronic device is used to provide computing and control capabilities. The memory of the electronic device includes a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system and a computer program. The internal memory provides an environment for the operation of the operating system and the computer program in the non-volatile storage medium. The communication interface of the electronic device is used to communicate with an external terminal in a wired or wireless manner. The wireless manner can be implemented through WIFI, a carrier network, near field communication (NFC), or other technologies. The display screen of the electronic device can be a liquid crystal display screen or an electronic ink display screen. The input device of the electronic device can be a touch layer covering the display screen, or a button, a trackball, or a touchpad provided on the housing of the electronic device, or an external keyboard, touchpad, or mouse, etc.

[0179] Those skilled in the art can understand that Figure 9 the structure shown in the figure is only a structural diagram of a part related to the technical solution of the present disclosure, and does not constitute a limitation on the electronic device to which the solution of this application is applied. The specific electronic device may include more or fewer components than those shown in the figure, or combine certain components, or have different component arrangements.

[0180] The fourth aspect of the present invention discloses a computer-readable storage medium. A computer program is stored on the computer-readable storage medium. When the computer program is executed by a processor, the steps in an inversion method of a base hybrid regularization model according to any one of the first aspects disclosed in the present invention are implemented.

[0181] Please note that the technical features of the above embodiments can be combined arbitrarily. For the sake of brevity of description, all possible combinations of the technical features in the above embodiments are not described. However, as long as the combinations of these technical features do not conflict, they should be considered as within the scope described in this specification. The above embodiments only represent several implementation manners of this application, and their descriptions are relatively specific and detailed, but they should not be construed as a limitation on the scope of the invention patent. It should be pointed out that for those of ordinary skill in the art, without departing from the concept of this application, several modifications and improvements can be made, and these all belong to the protection scope of this application. Therefore, the protection scope of the patent of this application should be subject to the appended claims.

[0182] The above is the preferred implementation manner of the present invention. It should be pointed out that for those of ordinary skill in the art, without departing from the principle of the present invention, several improvements and refinements can be made, and these improvements and refinements should also be regarded as within the protection scope of the present invention.

Claims

1. An inversion method based on a hybrid regularization model, characterized in that The method includes: Step S1, input seismic data and construct an objective function; Step S2, determine whether the number of iterations of the objective function reaches the maximum value Kmax and whether the error between two consecutive times reaches the minimum value; Step S3: Calculate the pseudo-gradient of the objective function according to the judgment result of Step 2 Step S4. Use L - BFGS to perform cyclic iteration to calculate the descent direction of the pseudo - gradient and map it to the same quadrant for correction; Step S5. Calculate a new inversion result through a linear search strategy according to the descent direction of the pseudo-gradient .

2. The inversion method based on the hybrid regularization model according to claim 1, characterized in that, In Step S1, when the seismic wave is vertically incident, the reflection coefficient of adjacent formation interfaces can be represented by the acoustic impedance of the formation. The formation acoustic impedance is continuous. Therefore, the reflection coefficient of the formation is expressed as: For a single seismic trace, the specific formula is: d = Gm + n where \(d\in R\) N represents single - trace seismic data, which is a vector with respect to the time series \(t\); \(G\) is the forward operator, composed of \(W*L\), \(W = R\) N×N is the convolution operator composed of source wavelets; \(L\) is the finite - difference operator, \(m\in R\) N represents the elastic parameter vector of a single - trace formation; \(n\in R\) N represents the random noise vector; Therefore, the objective function is:

3. The inversion method based on a base hybrid regularization model according to claim 2, characterized in that In Step S2, when the number of iterations of the objective function reaches the maximum value Kmax and the error between two consecutive times reaches a small value, stop the iteration; otherwise, continue the iteration.

4. The inversion method based on the base mixed regularization model according to claim 3, characterized in that In step S3, when the iteration in step S2 stops, calculate the pseudo-gradient of the objective function The calculation method is as follows: Among them, g i (m) is the gradient of the differentiable L2 norm term in the objective function, sign(*) represents the vector sign function, and ξ i represents the i-th element after vectorizing the L1 norm term, and β is the regularization parameter of the total variation term.

5. The inversion method based on the hybrid regularization model according to claim 4, characterized in that In Step S4, for solving large-scale non-linear optimization, calculate the inverse of the approximate Hessian matrix through L-BFGS, so as to obtain the descent direction of the pseudo-gradient, that is: where, H k is the inverse of the Hessian matrix of the objective function calculated by L - BFGS; is the pseudo - gradient of the objective function; k is the number of iterations; According to the descent direction of the pseudo-gradient, the vectors before and after can be updated to the same quadrant, and the orthogonal function is defined as: where sign(*) represents the vector sign function, and x i represents the formation parameter result of the first vector in the i-th layer; To ensure that the two variables fall in the same quadrant, the new descent direction of the pseudo-gradient is:

6. The inversion method based on a base hybrid regularization model according to claim 4, wherein In the step S5, the step size α of the solution result k The strategy adopted for linear search is as follows: Among them, When When The step size α k is the step size calculated by adopting a linear search strategy; represents the formation parameters of the i-th layer at the k-th iteration, represents one of the vectors of the orthogonal function, and sign(·) represents the sign of the vector.

7. The inversion method based on the hybrid regularization model according to claim 6, wherein In Step S5, when the number of iterations k reaches the maximum or the difference between the results of two consecutive iterations is within a certain range, end the algorithm and output the corresponding inversion result.

8. An inversion system for a base hybrid regularization model, characterized in that, The system includes: A first processing module, configured to input seismic data and construct an objective function; A second processing module, configured to determine whether the number of iterations of the objective function reaches the maximum value Kmax and whether the error between two consecutive times reaches a small value; A third processing module, configured to calculate a pseudo-gradient of the objective function according to the determination result of step 2 The fourth processing module is configured to calculate the pseudo-gradient by using L-BFGS cyclic iteration for the descending direction and map and correct the same quadrant; The fifth processing module is configured to calculate a new inversion result through a linear search strategy according to the descent direction of the pseudo-gradient .

9. An electronic device, characterized in that, The electronic device includes a memory and a processor. When the processor executes the computer program stored in the memory, it implements the steps in the inversion method of a base hybrid regularization model according to any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that, A computer program is stored on the computer-readable storage medium. When the computer program is executed by the processor, it implements the steps in the inversion method of a base hybrid regularization model according to any one of claims 1 to 7.

Images