CN110618453A

CN110618453A - Wave impedance inversion method based on improved damping least square method

Info

Publication number: CN110618453A
Application number: CN201910724007.0A
Authority: CN
Inventors: 吴朝容; 满建; 付小念; 段文燊
Original assignee: Chengdu Univeristy of Technology
Current assignee: Chengdu Univeristy of Technology
Priority date: 2019-08-07
Filing date: 2019-08-07
Publication date: 2019-12-27
Anticipated expiration: 2039-08-07
Also published as: CN110618453B

Abstract

A wave impedance inversion method based on an improved damping least square method belongs to the field of geophysical inversion, in particular relates to a wave impedance inversion technology in oil-gas geophysical exploration, and aims to provide an improved wave impedance inversion method which is used for establishing a wave impedance inversion target function in the least square sense, reducing the calculation cost of a Jacobian matrix in the traditional least square method, solving the problems of non-positive definite and singular of the Jacobian matrix, further optimizing the iteration mode of a damping coefficient, enabling the algorithm to have better stability and global convergence and enabling the inversion result to be more reliable, wherein the implementation process comprises the following steps: firstly, constructing an inversion target function according to a convolution model and a sparse constraint pulse inversion theory; establishing an initial model according to the logging information and the horizon information; setting algorithm and model parameters; putting the initial model into a target function for iteration; obtaining the optimal matching reflection coefficient under the comprehensive constraint conditions of three data volumes of model constraint, sparse constraint and seismic data constraint; solving the wave impedance by adopting a recursive method; and seventhly, adding high-low frequency compensation to obtain a broadband wave impedance value.

Description

Wave impedance inversion method based on improved damping least square method

Technical Field

The invention provides a wave impedance inversion based on an improved damping least square method, belongs to the field of geophysical inversion, researches a damping coefficient given mode under a confidence domain method, establishes a wave impedance inversion target function under the least square meaning, reduces the calculation cost of a Jacobian matrix in the traditional least square method, solves the problems of non-positive determination and singularity of the Jacobian matrix, optimizes an iteration mode of a damping coefficient, and enables an algorithm to have better stability and global convergence.

Background

Wave impedance inversion technology is continuously improved from the aspects of algorithm optimization, inversion strategies, modeling modes and the like to seek breakthrough from the beginning of development. The geophysical inverse problem is generally nonlinear, namely, no linear relation exists between observation data and model parameters, and the idea of linear inversion is to linearize a nonlinear problem by a nonlinear equation in inversion through the idea of parameter replacement, so that the aim of solving is fulfilled. The mature linear inversion algorithm such as steepest descent, conjugate gradient and least square obtains better results in the field of wave impedance inversion.

Least squares was first proposed by french mathematician Legendre (1806) in studies on the shape of planets and the gravity of spheres, and was first used to solve the linear fitting problem, and later, the method was further developed by scholars k.levenberg and d.marquardt in combination with a nonlinear two-multiplication solving problem to propose a damped least squares method, which is also called Levenberg-Marquardt algorithm, abbreviated as LM algorithm. The damping least square method is an improved form of the Gauss-Newton method, has the local convergence of the Gauss-Newton method and the global characteristic of the gradient descent method, and overcomes the requirement that the Gauss-Newton method has full rank for a Jacobian matrix. The damping least square method algorithm is mainly used for solving the optimization problem of the nonlinear multi-element objective function, and the current research on the damping least square method algorithm mainly focuses on the problems of convergence characteristics, stability and damping factor selection. For example, Osborne (1976) introduces an iterative weight factor to ensure the global convergence of the algorithm; more (1978) proposes a damping coefficient iteration rule based on a trust domain algorithm; chen Qi (1988) organically combines the damping least square method with the steepest descent method by adding a high-cut-off damping factor, thereby overcoming the high-damping factor bandThe disadvantages of the prior art; three-point search method is proposed by Chendhao et al (1994), which improves the calculation speed of the algorithm; yamashita and Fukuahima prove that when | | f (x) | I is locally error bounded, take mu_k＝||F(x_k)||²The algorithm has quadratic convergence; fan and Yuan (2012) demonstrate when the damping factor μ_k＝||F(x_k)||^δ，δ∈[1,2]The quadratic convergence of the time algorithm; lujin (2014) proposes an adaptive damped least squares method.

In the field of geophysics, the Ozong sea (1995) uses a damped least squares algorithm to calculate seismic source acceleration for an earthquake; kalachand et al (1995) applied the damped weighted least squares method to wide-angle seismic reflection time inversion; Ju-Wonoh et al (2013) apply damped least squares to a frequency domain full waveform inversion; YouzuoLin et al (2016) developed a damped least squares method based on Krift subspace, projected the original inverse problem to the Krift subspace in a dimensionality reduction manner, combined with parallel computing technology to obtain better results in a two-dimensional groundwater diffusion model; the Cainiello neural network algorithm and the damped least square method are combined and applied to porosity inversion by CyrilD, Boateng (2017) and the like, and convergence of the algorithm is improved.

The nonlinear equation in the geophysical inversion is often highly ill-conditioned or singular, and the damping coefficient introduced in the damping least square method can be understood to be regularized to some extent to solve the ill-conditioned condition of the equation set. The classical damped least square method has better robustness and faster convergence speed than the steepest descent method and the gauss-newton method, and is widely applied to the field of geophysical inversion all the time, however, whether the medium damping coefficient updating rule of the damped least square method is good or not directly influences the stability and the convergence of the algorithm. Testing for different damping coefficient values increases the computational cost for each iteration in the least squares sense. Therefore, it is urgent to develop a damped least squares algorithm that is not only computationally efficient, robust, and reduces the search cost. For most existing damped least square methods, the focus is mostly focused on solving the Jacobian matrix, orthogonal triangle (QR) decomposition or Singular Value Decomposition (SVD) is two direct solving methods, and the two methods have good stability for searching directions. However, the problem with both solutions is that they have too many parameters to solve, which is computationally expensive. The damping least square method based on the confidence domain method is adopted to convert the given problem of the damping coefficient into the subproblem of solving the confidence domain, the problems of local convergence, high Jacobian matrix calculation cost and unstable algorithm of the conventional linear inversion algorithm can be effectively solved, and the confidence domain method has the global convergence characteristic and high calculation efficiency.

Disclosure of Invention

Sparse constraint impulse inversion is mainly used for solving some problems existing in seismic record underdetermination, and the method generally assumes that the wave impedance of the underground geologic body and the layer reflection coefficient corresponding to the wave impedance are sparsely distributed, namely the wave impedance mainly comprises a stronger reflection coefficient sequence and a weaker Gaussian background reflection coefficient sequence. Based on this assumption, there are two inversion concepts.

The train of thought is that initial reflection coefficient is obtained by pulse deconvolution, maximum likelihood deconvolution or corresponding threshold processing technology, earthquake synthetic record is obtained by initial model and wavelet convolution, residual error of synthetic earthquake record and actual earthquake data in least square sense is obtained, and the number of reflection sequences is controlled and corrected by the residual error. The iterative cycle is adopted, an optimal reflection coefficient sequence can be finally obtained within an error range, and the aim of sparse pulse inversion is fulfilled.

Another idea is to synthesize the reflection coefficient constraint, seismic trace constraint and logging data constraint terms, and solve the overall minimum of the three constraint terms in the least square sense, and the objective function can be recorded as:

J＝∑|r_i|^p+λ^q∑(d_i-s_i)^q+α²∑(t_i-z_i)²→min

wherein: r is_i、λ、d_iRespectively representing a reflection coefficient, a sparse constraint coefficient and seismic data; alpha, s_i、t_i、z_iWeight coefficient and number of synthetic channels respectively representing data matchingAccording to the impedance constraint trend and the wave impedance constraint range; p, q represent L modulo factors, and p is usually 1 and q is 2. In the formula (1), Σ | r_i|^pRepresenting the sum of absolute values of the reflection coefficients, λ^q∑(d_i-s_i)^qRepresenting the difference, alpha, between the seismic data and the synthetic record²∑(t_i-z_i)²Is the sum of the squares of the differences in wave impedance trends.

In the calculation process, an initial model is selected firstly, the model is generally generated by a small number of pulses, and the model is modified through residual error control until the value of an objective function reaches the minimum. And then the control model gradually increases the number of pulses without changing, the previous iterative calculation is continued, the calculation is stopped until the inversion result is stable, and the inversion result is the final result.

The above solving process can also be expressed mathematically as solving an optimization problem with constraints, the expression is as follows:

in the formula: r is_iDenotes the value of the reflection coefficient, ZL_i,ZU_iRespectively representing the upper and lower limits of the inversion wave impedance constraint, d_i、λ、s_iRespectively expressed as actual seismic trace data, sparse constraint weight coefficients and convolution synthesized seismic traces.

The matching degree of the verification data and the reflection coefficient or the output impedance is adjusted by a sparse constraint weight coefficient lambda, and in addition, factors control the inversion effect, such as seismic wavelets, the constraint condition of the impedance and the weight coefficient alpha.

And solving according to the objective function to obtain an optimal reflection coefficient, and then obtaining the relative wave impedance by a trace integration or recursion method.

The classical damping least square method comprises the following implementation steps:

step 1: given an initial point x⁽¹⁾Growth factor beta > 1, initial parameter alpha₁> 0, allowable error ε > 0, calculate F (x)⁽¹⁾) α is set to α₁，k＝1。

Step 2: setting alpha as alpha/beta. Computing

And 3, step 3: solution equation

And 3, step 3: to find the direction d^(k)Let us order

x^(k+1)＝x^(k)+d^(k)

And 4, step 4: calculating F (x)^(k+1)) If F (x)^(k+1))＜F(x^(k)) And (6) turning to the step (6), otherwise, performing the step 5.

And 6, step 6: if it isStopping the calculation to obtain the solutionOtherwise, setting alpha to be beta alpha, and turning to the step (3).

And 7, step 7: if it isStopping the calculation to obtain a solutionOtherwise, k: ═ k +1 is set, and the step 2 is returned.

The strategy used for both the steepest descent method and the conjugate gradient method is to give a point x^(k)Then, a search direction d is defined^(k)From x again^(k)Starting edge d^(k)And performing one-dimensional search. For all inverse problems, all local optimization algorithms have the condition that the objective function falls into a local extremum, the search is finished, and the finally obtained solution is a local optimal solution instead of a global optimal solution.

The confidence domain method finds a new one aiming at the common fault of the linear optimization methodThe new solution, at the initial iteration point x^(k)Given this, a range of variation is determined, usually taken as x^(k)A sphere region as a center, called confidence region, the radius r of the sphere region, called confidence interval, the quadratic approximation formula of the objective function is optimized in the region, and the subsequent point x is obtained according to a certain mode^(k+1). If the precision does not meet the requirement, x is defined again^(k+1)And optimizing a new quadratic approximation formula for the objective function in the confidence domain of the center of the spherical domain until the precision requirement is met.

The unconstrained optimization problem is written as:

min f(x),x∈Rⁿ

for f (x) at the initial iteration point x^(k)And (3) second-order Taylor expansion:

wherein o ((x-x)^(k))²) A Pianogo remainder term, which is a high order infinitesimal representation, can generally be ignored during the solution process, and represents a Hamiltonian.

The above equation, with the higher order infinity term truncated, is approximately expressed as:

let d be x-x^(k)And (3) solving a quadratic function model:

limiting the value of d to x^(k)Nearby useApproximation f (x)^(k)+ d), making | | d | | | | less than or equal to r_k，r_kIs a given constant (confidence domain radius). The minimization problem of the function f (x) is attributed to the following sub-problems:

w is noted as:

then find d^(k)The requirements for being an optimal solution are:

let us²f(x^(k)) + wI is reversible, resulting from the above formula:

||d^(k)||＝||(▽²f(x^(k))+wI)^-1▽f(x^(k))||

from the above conditions, the solution d of the above formula^(k)And radius of confidence domain r_kIs related to the value of (A). If r is_kSufficiently large, the value of w is likely to be small, d^(k)Near newton direction, i.e.:

d^(k)≈-▽²f(x^(k))^-1▽f(x^(k))

if r_k→ 0, then | | d^(k)| → 0, w → + ∞, at which time:

i.e. d^(k)Approaching the steepest descent direction. Radius r of confidence domain_kWhen increasing gradually, d^(k)Continuously changing between the Newton direction and the steepest descent direction.

Obtain the optimal solution d^(k)Then, whether the approximate solution of the original unconstrained optimization problem can be viewed as point x^(k)+d^(k)Also according to the useApproaches f (x) success or not. Whether the success is successful or not can be determined according to the actual reduction and the pre-reduction of the function valueBy measuring the ratio of the decrease, i.e.

If ρ_kToo small a value, the successor point still takes x^(k)(ii) a If ρ_kIf the value is relatively large, the approximation is considered to be successful, and then x is made^(k+1)＝x^(k)+d^(k)。

Drawings

TABLE 1 three optimization algorithm parameter settings

Table 1 sets up the parameters for the three optimization algorithms. By using the parameters in table 1, the inversion of the relative wave impedance of the three optimization methods is performed under the conditions of unchanged wavelets and no noise, and the inversion results are shown in fig. 1, 2 and 3.

Fig. 1 is a wave impedance inversion result of a classical damped least square method, fig. 2 is a wave impedance inversion result of an adaptive damped least square method, and fig. 3 is a wave impedance inversion result of a confidence domain method. And analyzing the inversion accuracy of the three algorithms and the seismic channel residual error distribution condition under the condition that control parameters such as an error interval, iteration times and the like are not changed. The errors of the three inversion methods are respectively 1.3%, 0.9% and 0.8%, the inversion errors of the three inversion methods are not large, but the trust domain method is superior to the former two methods in the optimization strategy and has the characteristic of global convergence.

FIG. 4 is a normalized seismic record. The test of an inversion algorithm is carried out by utilizing a theoretical Marmousi model, the test model contains 1237 channels in total, the sampling time is 255ms, the sampling interval is 2ms, the spectrum analysis is carried out on the original seismic data, the recording frequency of the original seismic is 0-100HZ, the dominant frequency is about 30HZ, and according to the theory of seismic wave impedance inversion hypothesis, the zero-phase Rake wavelet with the inverted wavelet dominant frequency of 30HZ is assumed.

Actually, the target function of the following formula is adopted for inversion, and the seismic record shown in fig. 4 is obtained after normalization processing is carried out on the result after inversion.

J＝∑|r_i|^p+λ^q∑(d_i-s_i)^q+α²∑(t_i-z_i)²→min

Fig. 5 is the normalized initial model. Since there is no horizon and logging information, the low frequency model in FIG. 5 is obtained by smoothing the original wave impedance profile. After algorithm optimization parameters such as damping coefficients, iteration termination parameters and sparse control coefficients are selected, the influence of model constraint parameters on inversion results is researched under the condition that the parameters are not changed, and part of inversion parameters are shown in table 2.

TABLE 2 related inversion parameter settings

Parameter(s)	Value taking
		Sparse constraint coefficients	0.01
Initial trust domain radius	0.1
		Wavelet	Zero-phase Rake wavelet with 50HZ main frequency

Fig. 6 to 11 are relative wave impedance profiles obtained by gradually changing the model constraint coefficient α and inverting the model constraint coefficient α by using an improved damped least squares method under the premise that the parameters in table 2 are not changed.

Fig. 6 is a relative wave impedance profile obtained by inversion when α is 0, fig. 7 is a relative wave impedance profile obtained by inversion when α is 0.001, fig. 8 is a relative wave impedance profile obtained by inversion when α is 0.05, fig. 9 is a relative wave impedance profile obtained by inversion when α is 0.1, fig. 10 is a relative wave impedance profile obtained by inversion when α is 1, and fig. 11 is a relative wave impedance profile obtained by inversion when α is 10. As can be seen from fig. 6, 7, 8, 9, 10, and 11, when the model constraint coefficient α is 0, the inverted relative wave impedance profile has a high longitudinal resolution, but a poor transverse continuity; when alpha is 0.001, even if the weight of the model constraint term is smaller, the transverse continuity of the reflection coefficient profile and the relative wave impedance profile obtained by inversion is obviously improved, and the important function of the model constraint term on the whole inversion result is satisfied; gradually increasing the value of the model constraint coefficient, as shown in fig. 8, 9, 10 and 11, when α is sequentially increased from 0.05 to 0.1, 1 and 10, the lateral continuity of the inversion section is gradually increased; when the weight of the model constraint item is too large, the corresponding seismic data constraint action is relatively weak, the inversion result is excessively modeled, and a high-resolution wave impedance inversion section cannot be obtained.

Detailed Description

The trusted domain method comprises the following implementation steps:

step 1: given confidence domain radius r₁Feasible point x⁽¹⁾Accuracy requirement_εAnd parameters

0＜μ＜η＜1 (general getting)) And k: ═ 1.

Step 2: by calculating f (x)^(k)),▽f(x^(k)) If | | f (x)^(k)) If | | < epsilon, stopping the calculation to obtain x^(k)(ii) a Otherwise, calculate ^²f(x^(k)).

And 3, step 3: solving sub-problems

Solving the optimal solution d of the subproblem^(k). Order:

and 4, step 4: if ρ_kLess than or equal to mu, let x^(k+1)＝x^(k)(ii) a If ρ_kMu, let x^(k+1)＝x^(k)+d^(k).

And 5, step 5: modifying r_kIf ρ_kLess than or equal to mu, makeIf μ < ρ_kEta, r_k+1＝r_k(ii) a If ρ_kGreater than or equal to eta, let r_k+1＝2r_k.

And 6, step 6: and (5) setting k to be k +1, and turning to the step (2).

Claims

1. A wave impedance inversion method based on an improved damped least square method is characterized in that a set of brand-new and effective wave impedance inversion method is formed, and comprises the following steps:

(1) constructing a wave impedance inversion target function according to a convolution model and a sparse constraint pulse theory:

J＝∑|r_i|^p+λ^q∑(d_i-s_i)^q+α²∑(t_i-z_i)²→min \*MERGEFORMAT(1

in the formula, r_i、λ、d_iRespectively representing a reflection coefficient, a sparse constraint coefficient and seismic data; alpha, s_i、t_i、z_iRespectively representing the weight coefficient of data matching, synthetic channel data, impedance constraint trend and wave impedance constraint range; p, q represent modulo factors, usually taken as L. In the formula (1), Σ | r_i|^pRepresenting the sum of absolute values of the reflection coefficients, λ^q∑(d_i-s_i)^qRepresenting the difference, alpha, between the seismic data and the synthetic record²∑(t_i-z_i)²Is the sum of the squares of the differences in wave impedance trends.

(2) An initial model is selected, which is generally generated by a smaller number of pulses, and the model is modified by residual control until the value of the objective function is minimized. And then the control model gradually increases the number of pulses without changing, the previous iterative calculation is continued, the calculation is stopped until the inversion result is stable, and the inversion result is the final result.

(3) On the basis of establishing an optimization objective function and an initial model, an iterative optimization rule based on a trust domain method is established, and the specific implementation mode is as follows:

step 1: given confidence domain radius r₁Feasible point x⁽¹⁾Accuracy requirement_εAnd the parameter 0 < mu < eta < 1 (generally taken)) And k: ═ 1.

Step 2: by calculating f (x)^(k)),If it isStopping the calculation to obtain x^(k)(ii) a Otherwise, calculating

And 3, step 3: solving sub-problems

s.t.||d||≤r_k,

Solving the optimal solution d of the subproblem^(k). Order:

And 6, step 6: and (5) setting k to be k +1, and turning to the step (2).

(4) And obtaining the reflection coefficient under the least square meaning according to the third step.

(5) And converting the reflection coefficient in the fourth step into relative wave impedance by adopting a recursive method or a trace integration method.

(6) And adding low-frequency and high-frequency models into the relative wave impedance in the fifth step to obtain the broadband wave impedance.

2. The wave impedance inversion method based on the improved damped least squares method according to claim 1, wherein the wave impedance inversion method comprises the following steps: compared with the objective function established by the fluctuation equation based on the fluctuation theory, the objective function established on the basis of the convolution model has the characteristics of strong anti-noise performance and stable algorithm.

3. The wave impedance inversion method based on the improved damped least squares method according to claim 1, wherein the wave impedance inversion method comprises the following steps: and the sparse constraint pulse inversion theory which is widely applied in inversion is combined, the inversion result is faithful to seismic data, and the multi-solution is good.

4. The wave impedance inversion method based on the improved damped least squares method according to claim 1, wherein the wave impedance inversion method comprises the following steps: the damping least square method based on the confidence domain algorithm is applied to the post-stack wave impedance inversion method for the first time, compared with the conventional linear inversion algorithm such as the steepest descent method, the conjugate gradient method, the Gauss Newton method and the like, the confidence domain method breaks through the conventional optimization strategy, introduces the concept of the confidence domain, and effectively improves the robustness, the calculation efficiency and the global convergence of the algorithm.