CN114966831A - Viscoacoustic full waveform inversion method based on velocity-attenuation decoupling - Google Patents

Abstract

The invention belongs to the technical field of seismic exploration, and particularly relates to a visco-acoustic full waveform inversion method based on velocity-attenuation decoupling, which comprises the following steps of: 1. setting seismic source, wave detection point and seismic record information, and reading observation records; 2. solving a fractional order viscoacoustic equation according to the time forward direction to obtain simulation data and storing a part of wave fields; 3. and if the maximum iteration times are not reached, repeating the step two until the maximum iteration times are reached and obtaining the reconstruction speed and Q. The speed and Q parameters with higher resolution can be recovered; the decoupling technology based on TRN can effectively relieve the coupling effect in multi-parameter simultaneous inversion and improve the inversion precision; the truncation condition in the TRN algorithm can automatically judge whether to obtain an accurate descending direction, so that over-updating or under-updating is avoided, and the algorithm complexity is reduced; the method can effectively inhibit crosstalk, can better recover the speed and attenuation model, and can meet the requirements of the oil and gas industry on the parameters of the underground medium.

Description

Viscoacoustic full waveform inversion method based on velocity-attenuation decoupling

Technical Field

The invention belongs to the technical field of seismic exploration, and particularly relates to a velocity-attenuation decoupling-based viscoacoustic full waveform inversion method.

Background

Oil and gas resources are black gold in the industry and are also important strategic resources of the country. At present, deep layer oil gas has become an important development field of global petroleum industry due to continuous deep exploration and development. However, deep seismic waves have long propagation paths and strong signal attenuation. Thus. Deep oil and gas reservoirs are difficult to accurately identify, and the exploration and development risks are large. Therefore, obtaining accurate seismic wave propagation velocity and quality factor Q of the stratum is an important prerequisite for accurately identifying deep oil and gas reservoirs. The existing speed and Q modeling technology cannot obtain high-precision speed and Q fields at the same time. Therefore, there is a need for a method of reconstructing subsurface parameters that can achieve both high precision velocity and Q.

The full waveform inversion method (FWI) is a method of reconstructing the parameters of the subsurface medium by matching all the valid data in the received records, and is theoretically the most accurate modeling method at present. In contrast to acoustic FWI, the sticky acoustic full waveform inversion (QFWI) can compensate for energy attenuation during seismic wave propagation and can simultaneously invert the velocity and Q-field of complex media. Chen et al (2017) studied a full waveform inversion method based on fractional order ordinary Q viscous acoustic wave equation and proposed the feasibility of velocity-attenuation two-parameter inversion. Scott et al (2018) apply a multi-scale inversion strategy in the full waveform inversion of the sticky sound waves, and effectively reduce the problem of crosstalk between speed and Q by using a Gaussian Newton method. Xue et al (2018) propose a gradient precondition depending on Q, accelerate the convergence speed of time domain viscoacoustic full-waveform inversion, compensate the loss of Q gradient by inverting the sign of amplitude loss, and reduce the sensitivity of Q attenuation. However, velocity-attenuation bi-parametric inversion has not been an effective strategy. This is because although in theory QFWI can achieve both velocity and Q fields with high precision. But since the seismic response due to velocity perturbations and Q perturbations are similar, this strong coupling effect causes the velocity update and Q update to interfere with each other during the inversion process, which is also commonly referred to as a cross-talk problem between parameters. Currently, most FWIs update parameters using local optimization methods based on gradient classes (e.g. conjugate gradient method). Although the calculation amount is relatively low, the number of iterations required for achieving convergence is very high, and the problem of crosstalk between different parameters is difficult to alleviate. In addition, the inversion accuracy is seriously influenced by the strong coupling effect of the velocity-attenuation double parameters.

Meanwhile, the problems of data band limitation, uneven illumination, multi-parameter coupling and the like can be effectively relieved by utilizing the gradient and Heisen information for inversion. However, directly solving the hessian matrix's inverse is very computationally intensive (e.g., if one wants to solve exactly a two-dimensional model of size nx × nz, one needs to do nx × nz ^nx×nz Paradoxical), current computing hardware is still unable to meet this requirement. Therefore, there is a need for an efficient optimization scheme to obtain low cost, high precision, approximately hessian operators. Operto proposes a Gaussian-Newton method to approximate the Hessian operator, but the method has high requirements on calculation amount and storage. Finite memory storage (BFGS) constructs the approximate hessian operator using gradients in an iterative process. However, this method consumes a huge amount of memory, and is not suitable for processing large-scale data. In order to save memory consumption, the L-BFGS estimates approximate Hessian operators through the last few models and gradient information, however, the matrix obtained by the method is positive definite, and the correct Hessian operator is not positive definite. For highly ill-conditioned FWI, this positive approximation may lead to undesirable inversion results. Unlike the above method, the truncated newton method (TRN) does not require explicit solution of the hessian operator, but only requires iterative updating of the inaccurate (truncated) solution associated with the newton equation by a small number of conjugate gradients at each iteration of FWI. The approximate hessian operator can be constructed at a lower cost.

Disclosure of Invention

In order to reconstruct high-precision speed and Q parameters, the invention provides a viscoacoustic full waveform inversion method based on speed-attenuation decoupling. Aiming at a fractional order visco-acoustic equation, a dual-parameter QFWI for simultaneously inverting speed and Q is provided; the truncated Newton method is introduced, the Hessian operator is effectively approximated, the problems of speed and Q strong coupling effect are relieved, and the dual-parameter inversion precision is improved; a judgment condition is introduced into a TRN algorithm to avoid insufficient or excessive updating of the descending direction; the method can effectively inhibit crosstalk, can better recover the speed and attenuation model, and can meet the requirements of the oil and gas industry on the parameters of the underground medium.

The technical scheme adopted by the invention is as follows: a visco-acoustic full waveform inversion method based on velocity-attenuation decoupling comprises the following steps:

step one, setting seismic source, wave detection point and seismic record information, and reading observation data;

solving a fractional order viscoacoustic equation according to the time forward direction to obtain simulation data (formula 2), and storing a part of wave fields; acquiring an accompanying seismic source by using a formula 9, and reversely acquiring an accompanying wave field according to time by using a formula 7; reversely reconstructing a forward wave field by utilizing a part of wave fields stored in the flow 2 and solving a formula 5 and a formula 6; solving the speed and the Q gradient by using a formula 4; the descending direction is solved by using an algorithm 1; updating the model parameters using equation 10;

the gradient of the objective function is found, and for QFWI, the least squares objective function can be expressed as:

in the formula, m represents a model parameter in a space X; d and p represent observed and simulated data wavefields, x, respectively _r Representing the geophone position;

the simulated wavefield p may be calculated by seismic forward modeling, and the viscoacoustic equation with fractional laplace operator may be expressed as:

Ap＝f (2)；

where f denotes the seismic source, c denotes the medium velocity,

ω ₀ for reference to angular frequency, operator A is a Partial Differential Equation (PDEs) that describes the propagation of waves in the subsurface, the equationThe speed and Q in the process can both show expression, and the included parameters are fewer, so that the calculation of the gradient is facilitated. Using the lagrange multiplier method, the gradient G based on the objective function (1) can be expressed as:

where <, > represents the inner product, μ represents the adjoint wavefield, and the second term in the inner product can be represented as:

the adjoint wavefield of the first term in the inner product can be solved by:

in the formula (I), the compound is shown in the specification,

indicates the accompaniment of f ^* The representation is accompanied by a seismic source,

the method represents an adjoint operator which determines the scattering mode of model parameters, different wave equation systems and different characterization parameters, the scattering modes are different, and the adjoint operator in the method can be expressed as follows:

from equation 7, the expression form accompanying the seismic source is related to the form of the objective function, and is expressed as:

f ^* ＝p-d (9)

since QFWI can be considered as a PDE constraint optimization problem, the update model parameter m can be expressed as:

m _k+1 ＝m _k +a _k d _k (10)

in the formula (d) _k And alpha _k Respectively representing the descending direction and descending step length at the k-th iteration, the step length can be estimated by linear search, and the model updating quantity d _k The solution can be performed by a gradient-like method (equation 11) or a newton-like method (equation 12):

d _k ＝-P(m _k )G(m _k ) (11)

d _k ＝-H ^-1 (m _k )G(m _k ) (12)

wherein P represents a precondition, G represents a gradient, and H represents a Hessian matrix;

the gradient-like method (equation 11) that ignores the hessian matrix effect, although the amount of computation is relatively low, the number of iterations required to achieve convergence is very high; although the convergence speed can be accelerated by using a precondition method such as L-BFGS (L-BFGS), the decoupling effect of the precondition on multi-parameter crosstalk is limited; the Newton method (formula 12) can relieve the problems of data band limitation, uneven illumination, multi-parameter coupling and the like, and can improve the convergence of inversion; however, most of the methods need to obtain the inverse of the hessian matrix, the calculation amount is huge, and the current calculation hardware cannot meet the requirement;

therefore, the problem of crosstalk of simultaneous inversion of multiple parameters in QWI is solved; the second key technology of the method is to approximately solve the Newton descending direction by using TRN; the method can avoid directly solving the inverse of a Hessian matrix under the condition of no matrix, and a flow chart for solving the descending direction by utilizing TRN is shown as an algorithm 1;

the algorithm is as follows:

the algorithm gradually obtains the Newton's descending direction in an internal conjugate gradient cycle by introducing the product of the Hessian operator and a vector lambda (Hessian vector product, H lambda); with the second order adjoint state method, the hessian vector product H λ can be expressed as:

in the formula, the ratio of omega,

respectively representing a calculation domain and a model parameter, H (m) lambda represents a Hessian vector product of m parameters, and in order to avoid over solving a Newton descending direction while ensuring algorithm convergence, a truncation condition is set:

||H(m)λ+G(m)||＞0.9||G(m)|| (16)

furthermore, when a curvature related to the operator negative eigenvalue is encountered, i.e. α ₂ In the case of ≦ 0, the inner conjugate gradient cycle must be stopped, the descent direction being the result of the last inner conjugate gradient. The truncation condition ensures that the model updating amount of the inner conjugate gradient cyclic calculation is always ensured to be the descending direction of the objective function;

and step three, if the maximum iteration times are not reached, repeating the step two until the maximum iteration times are reached and obtaining the reconstruction speed and Q.

The invention has the beneficial effects that: a visco-acoustic full waveform inversion method based on velocity-attenuation decoupling is provided. Different from the traditional sound wave FWI, the algorithm of the method uses a constant fractional order wave equation as a wave equation system, and can simultaneously reconstruct speed and an attenuation parameter (Q); in order to avoid the crosstalk problem of multi-parameter simultaneous inversion, a local optimization strategy based on TRN is designed, the optimization strategy can avoid directly solving a hessian operator, the hessian vector product is efficiently calculated in a matrix-free mode by using a second-order adjoint state method, the descending direction of the model parameters is obtained by using a Newton-conjugate gradient frame, and an effective conjugate gradient termination condition is introduced in order to avoid the over-updating of the descending direction; the speed and Q parameters with higher resolution can be recovered; the decoupling technology based on TRN can effectively relieve the coupling effect in multi-parameter simultaneous inversion and improve the inversion accuracy; the truncation condition in the TRN algorithm can automatically judge whether the accurate descending direction is obtained, so that over-updating or under-updating is avoided, and the algorithm complexity is reduced; the method can effectively inhibit crosstalk, can better recover the speed and attenuation model, and can meet the requirements of the oil and gas industry on the parameters of the underground medium.

Drawings

FIG. 1 is a flow chart of a visco-acoustic full waveform inversion method based on velocity-attenuation decoupling;

FIG. 2 is a differential method velocity gradient plot of gradient validation;

FIG. 3 is a companion velocity gradient map for gradient validation;

FIG. 4 is a velocity gradient residual map of gradient validation;

FIG. 5 is a differential method Q gradient of gradient validation;

FIG. 6 is a companion Q gradient map of gradient validation;

FIG. 7 is the Q gradient residual of the gradient validation;

FIG. 8 is a differential method velocity plot for Hessian vector product validation;

FIG. 9 is a adjoint velocity plot for Hessian vector product validation;

FIG. 10 is a velocity residual plot for Hessian vector product validation;

FIG. 11 is a differential method Q plot of Hessian vector product validation;

FIG. 12 is a graph of the adjoint method Q of Hessian vector product validation;

FIG. 13 is a graph of the Q residual for Hessian vector product validation;

FIG. 14 is a diagram of a true velocity model;

FIG. 15 is a true Q model diagram;

FIG. 16 is an initial velocity model diagram;

FIG. 17 is an initial Q model diagram;

FIG. 18 is a gradient method based velocity pattern;

FIG. 19 is based on gradient method Q patterns;

FIG. 20 is a TRN decoupled speed pattern;

figure 21 is a TRN decoupled Q pattern;

FIG. 22 is a graph of the results of TRN-based velocity inversion;

FIG. 23 is a graph of the results of TRN-based Q-inversion;

FIG. 24 is a graph of the results of gradient-based velocity inversion;

fig. 25 is a graph of the results of gradient-based Q inversion.

Detailed Description

Example one

Referring to fig. 1, a viscoacoustic full waveform inversion method based on velocity-attenuation decoupling includes the following steps:

step one, setting seismic source, wave detection point and seismic record information, and reading observation data;

solving a fractional order viscoacoustic equation according to the time forward direction to obtain simulation data (formula 2), and storing a part of wave fields; acquiring an accompanying seismic source by using a formula 9, and reversely acquiring an accompanying wave field according to time by using a formula 7; reversely reconstructing a forward wave field by utilizing a part of wave fields stored in the flow 2 and solving a formula 5 and a formula 6; solving the speed and the Q gradient by using a formula 4; the descending direction is solved by using an algorithm 1; updating the model parameters using equation 10;

the gradient of the objective function is found, and for QFWI, the least squares objective function can be expressed as:

in the formula, m represents a model parameter in a space X; d and p represent observed and simulated data wavefields respectively,x _r representing the geophone position;

the simulated wavefield p may be calculated by seismic forward modeling, and the viscoacoustic equation with fractional laplace operator may be expressed as:

Ap＝f (2)；

where f denotes the seismic source, c denotes the medium velocity,

ω ₀ for reference to angular frequency, operator a is a Partial Differential Equation (PDEs) describing the propagation of waves in the ground, where both velocity and Q can be expressed, and the included parameters are few, which facilitates the calculation of the gradient, and the gradient G based on the objective function (1) can be expressed as:

where <, > represents the inner product, μ represents the adjoint wavefield, and the second term in the inner product can be represented as:

the adjoint wavefield for the first term in the inner product may be solved by:

in the formula (I), the compound is shown in the specification,

indicates the accompaniment of f ^* The representation is accompanied by a seismic source,

the method represents an adjoint operator which determines the scattering mode of model parameters, different wave equation systems and different characterization parameters, wherein the scattering mode is different, and the adjoint operator in the method can be expressed as follows:

from equation 7, the expression form accompanying the seismic source is related to the form of the objective function, and is expressed as:

f ^* ＝p-d (9)

since QFWI can be considered as a PDE constraint optimization problem, the update model parameter m can be expressed as:

m _k+1 ＝m _k +a _k d _k (10)

in the formula (d) _k And alpha _k Respectively representing the descending direction and descending step length at the k-th iteration, the step length can be estimated by linear search, and the model updating quantity d _k The solution can be performed by a gradient-like method (equation 11) or a newton-like method (equation 12):

d _k ＝-P(m _k )G(m _k ) (11)

d _k ＝-H ^-1 (m _k )G(m _k ) (12)

wherein P represents a precondition, G represents a gradient, and H represents a Hessian matrix;

the gradient-like method (equation 11) that ignores the hessian matrix effect, although the amount of computation is relatively low, the number of iterations required to achieve convergence is very high; although the convergence speed can be accelerated by using a precondition method such as L-BFGS (L-BFGS), the decoupling effect of the precondition on multi-parameter crosstalk is limited; the Newton method (formula 12) can relieve the problems of data band limitation, uneven illumination, multi-parameter coupling and the like, and can improve the convergence of inversion; however, most of the methods need to obtain the inverse of the hessian matrix, the calculation amount is huge, and the current calculation hardware cannot meet the requirement;

therefore, the problem of crosstalk of simultaneous inversion of multiple parameters in QWI is solved; the second key technology of the method is to approximately solve the Newton descending direction by using TRN; the method can avoid directly solving the inverse of a Hessian matrix under the condition of no matrix, and a flow chart for solving the descending direction by using TRN is shown as an algorithm 1;

the algorithm is as follows:

the algorithm gradually obtains the Newton's descending direction in an internal conjugate gradient cycle by introducing the product (Hessian vector product, H lambda) of a Hessian operator and a vector lambda; with the second order adjoint state method, the hessian vector product H λ can be expressed as:

in the formula, the reaction is divided into omega,

respectively representing a calculation domain and model parameters, wherein H (m) lambda represents a Hessian vector product of m parameters, and in order to avoid over solving the Newton descent direction while ensuring the convergence of the algorithm, a truncation condition is set:

||H(m)λ+G(m)||>0.9||G(m)|| (16)

in addition to this, the present invention is,when a curvature associated with the operator negative eigenvalue is encountered, i.e. alpha ₂ In the case of ≦ 0, the inter-conjugate gradient cycle must be stopped, with the descent direction being the result of the last inter-conjugate gradient. The truncation condition ensures that the model updating amount of the inner conjugate gradient cyclic calculation is always ensured to be the descending direction of the objective function;

and step three, if the maximum iteration times are not reached, repeating the step two until the maximum iteration times are reached and obtaining the reconstruction speed and Q.

Example two

In order to verify the correctness of the gradient and Hessian operator calculation in the method, a uniform velocity model is designed for verification, and the sizes in the horizontal direction and the vertical direction and the grid spacing are respectively 300m and 10 m. Ricker wavelets with a main frequency of 15Hz are used as excitation sound sources, the time sampling rate is 1ms, the real speed and the initial speed are 2050m/s and 2000m/s, and the real speed and the initial Q are respectively 80 and 100. For small perturbations of the model parameters, the gradient of the objective function to the model parameters can be calculated according to the definition of the gradient (difference method, equation 17):

in the formula (I), the compound is shown in the specification,

since the gradient of the ith point obtained by the definition method is expressed, the gradient of all the points can be obtained by performing the wave field simulation point by the definition method, and the gradient of all the points can be calculated by only 2 wave field simulations according to the state method, wherein the gradient is Δ m _i Representing a perturbation parameter.

In the experiment, the disturbance speed and the disturbance Q are respectively 1m/s and 20, and in order to observe a first Fresnel zone between the source and the detector, the source and the detector are respectively placed at models (x is 50m, y is 50m) and (x is 250m, y is 250 m). As shown in fig. 2-7, all are shown at the same scale. The gradients generated by the two methods are almost consistent, and the velocity gradient is abnormal only at the position of the seismic source, so the residual error of the velocity gradient can be ignored. Since the perturbation of Q is relatively large, there is a small amount of residual error in the Q gradient. But if the perturbation of Q continues to be reduced, the gradient energy of Q will be drowned in the numerical noise. However, as shown in fig. 5 and 6, the similarity is sufficient to verify the validity of our companion state method for gradient calculations. From the definition of the hessian matrix, it is known that the hessian matrix is equal to the derivative of the gradient:

in the formula, H _i，j (m) represents the hessian matrix element at the (i, j) position in the model, and when the vector λ in H (m) λ is (0,.. 1.. 0), H (m) λ ≈ H _i，j (m)。

Validation of the hessian vector product as shown in fig. 8-13, the perturbation point was placed in the center of the model. As can be seen from fig. 8 to 10, the result of the hessian vector product associated with the state method is highly consistent with the difference method, and only two outliers exist at the positions of the seismic source and the wave builder; and the Hessian vector product of Q has a disturbance anomaly at the disturbance point. But also sufficient to demonstrate the accuracy of the hessian vector product.

EXAMPLE III

In multi-parameter simultaneous inversion, the problem to be solved first is the crosstalk problem between different parameters. We verify the decoupling effect of Hessian vector product in multi-parameter inversion through an anomalous body model, the size of the model is 150 x 150 grid points, and the space intervals are all 10 m. The velocity model and the Q model are shown in FIGS. 14-17. FIGS. 14-17 represent the true velocity model, the true Q model, the initial velocity model, and the initial Q model, respectively. Ricker wavelets with a main frequency of 15Hz are used as excitation sound sources, and the time sampling rate is 1 ms. The 15 cannons are uniformly distributed on the ground surface, and the detectors are respectively distributed on the periphery of the model. The falling direction of the first iteration is shown in fig. 18-21, fig. 18 and 19 are based on the gradient class method, and fig. 20 and 21 are based on the TRN method. As can be seen from fig. 18 and 19, the gradient-based method has different degrees of crosstalk between different parameters in the descending direction. Although the anomaly homing in the speed-down direction is relatively accurate; but in the attenuation fall direction the velocity anomaly is erroneously mapped into the attenuation direction, even with energy far exceeding the attenuation anomaly. If the direction is directly used for parameter updating, the Q parameter may not be correctly inverted. In fig. 20 and 21, it can be observed that the energy of the abnormal body with different parameters is more focused, the problem of crosstalk between different parameters is effectively alleviated, and the attenuation abnormal body is accurately positioned. Inversion results obtained for different descent directions (10 iterations) are shown in fig. 22-25, where fig. 22-25 represent TRN-based velocity inversion results, respectively; q inversion result based on TRN; a gradient-based velocity inversion result; the result of the gradient-based Q inversion. It can be seen from an examination of fig. 22-25 that the velocity parameters reconstructed by the gradient-based method include a weak attenuation perturbation (fig. 24), while the attenuation parameters in fig. 25 include a large number of velocity "footprints" that severely interfere with the correct reconstruction of the attenuation parameters. Whereas in this method both velocity and attenuation anomalies are well reconstructed and there is no apparent cross talk problem (fig. 22 and 23). Therefore, the method can effectively suppress crosstalk among different parameters in QWI and can better recover the speed model and the attenuation model.

Claims (1)

1. A viscoacoustic full waveform inversion method based on velocity-attenuation decoupling is characterized in that: the method for inverting the viscoacoustic full waveform comprises the following steps:

step one, setting seismic source, wave detection point and seismic record information, and reading observation data;

solving a fractional order viscoacoustic equation according to the time forward direction to obtain analog data; solving an adjoint wave field; constructing a speed and a Q gradient; solving a descending direction and updating model parameters;

the gradient of the objective function is found, and for QFWI, the least squares objective function can be expressed as:

in the formula, m represents a model parameter in a space X; d and p represent observed and simulated data wavefields, x, respectively _r Representing the position of the detector；

The simulated wavefield p may be calculated by seismic forward modeling, and the viscoacoustic equation with fractional laplace operator may be expressed as:

Ap＝f (2)

where f denotes the seismic source, c denotes the medium velocity,

ω ₀ for reference to angular frequency, operator a is a partial differential equation describing the propagation of a wave in the subsurface, and using the lagrange multiplier method, the gradient G based on the objective function (1) can be expressed as:

where <, > represents the inner product and μ represents the adjoint wavefield, the second term in the inner product can be represented as:

the adjoint wavefield of the first term in the inner product can be solved by:

in the formula (I), the compound is shown in the specification,

indicates the accompaniment of f ^* The representation is accompanied by a seismic source,

the adjoint operator is represented, and the adjoint operator in the method can be expressed as:

the satellite sources are represented as:

f ^* ＝p-d (9)

the update model parameter m can be expressed as:

m _k+1 ＝m _k +a _k d _k (10)

in the formula (d) _k And alpha _k Respectively representing a descending direction and a descending step length in the kth iteration, wherein the step length can be estimated through linear search, and the descending direction is solved by using TRN;

the algorithm is as follows:

gradually acquiring a Newton descent direction in an internal conjugate gradient cycle by introducing a product of a Hessian operator and a vector lambda; with the second order adjoint state method, the hessian vector product H λ can be expressed as:

in the formula, the ratio of omega,

respectively representing a calculation domain and a model parameter, H (m) lambda represents a Hessian vector product of m parameters, and in order to avoid over solving a Newton descending direction while ensuring algorithm convergence, a truncation condition is set:

||H(m)λ+G(m)||>0.9||G(m)|| (16)；

and step three, if the maximum iteration times are not reached, repeating the step two until the maximum iteration times are reached and obtaining the reconstruction speed and Q.

Images