A viscoacoustic full waveform inversion method based on velocity-attenuation decoupling

Abstract

The present invention belongs to the field of seismic exploration technology, and specifically relates to a viscoacoustic full waveform inversion method based on velocity-attenuation decoupling, comprising the following steps: 1. Setting the source, detection point, seismic record information, and reading the observation record; 2. Solving the fractional-order viscoacoustic equation in the time forward direction, obtaining simulation data, and storing part of the wave field; 3. Repeating step 2 until the maximum number of iterations is reached and the reconstructed velocity and Q are obtained if the maximum number of iterations is not reached. The velocity and Q parameters with high resolution can be restored; the decoupling technology based on TRN can effectively alleviate the coupling effect in the simultaneous inversion of multiple parameters and improve the inversion accuracy; the truncation condition in the TRN algorithm can automatically determine whether the accurate descent direction is obtained, avoid over-update or insufficient update, and reduce the algorithm complexity; the method can effectively suppress crosstalk, and can better restore the velocity and attenuation model, and can meet the requirements of the oil and gas industry for underground medium parameters.

Description

Technical Field

[0001] The invention belongs to the technical field of seismic exploration, and in particular relates to a viscoacoustic full waveform inversion method based on velocity-attenuation decoupling. Background Art

[0002] Oil and gas resources are the "black gold" of the industry and an important strategic resource for the country. At present, due to the continuous deepening of exploration and development, deep oil and gas have become an important development field of the global oil industry. However, due to the long propagation path of deep seismic waves and strong signal attenuation. Therefore. Deep oil and gas reservoirs are difficult to accurately identify, and the risk of exploration and development is relatively high. Therefore, obtaining accurate seismic wave propagation velocity and the quality factor Q of the formation is an important prerequisite for accurately identifying deep oil and gas reservoirs. Existing velocity and Q modeling technologies cannot simultaneously obtain high-precision velocity and Q fields. Therefore, there is an urgent need for an underground parameter reconstruction method that can simultaneously obtain high-precision velocity and Q.

[0003] Full waveform inversion (FWI) is a method of reconstructing underground medium parameters by matching all valid data in the received records. In theory, it is the most accurate modeling method currently. Compared with acoustic FWI, viscoacoustic full waveform inversion (QFWI) can compensate for the energy attenuation during seismic wave propagation and can simultaneously invert the velocity and Q field of complex media. Chen et al. (2017) studied the full waveform inversion method based on the fractional-order constant Q viscoacoustic wave equation and proposed the feasibility of velocity-attenuation dual parameter inversion. Scott et al. (2018) applied a multi-scale inversion strategy in viscoacoustic full waveform inversion and used the Gauss-Newton method to effectively reduce the crosstalk problem between velocity and Q. Xue et al. (2018) proposed a Q-dependent gradient precondition to accelerate the convergence speed of viscoacoustic full waveform inversion in the time domain, and compensated for the loss of Q gradient by reversing the sign of amplitude loss, reducing the sensitivity of Q attenuation. However, there is no effective strategy for velocity-attenuation dual parameter inversion. This is because although in theory, QFWI can simultaneously obtain high-precision velocity field and Q field. However, since the seismic responses caused by velocity disturbances and Q disturbances are similar, this strong coupling effect causes the update of velocity and the update of Q to interfere with each other during the inversion process, which is often referred to as the crosstalk problem between parameters. At present, most FWI uses gradient-based local optimization methods to update parameters (such as the conjugate gradient method). Although this method has a relatively low computational complexity, the number of iterations required to reach convergence is very high, and it is difficult to alleviate the crosstalk problem between different parameters. In addition, the strong coupling of the velocity-attenuation dual parameters seriously affects its inversion accuracy.

[0004] Using gradient and Hessian information for inversion can effectively alleviate problems such as data band limitation, uneven illumination, and multi-parameter coupling. However, the amount of calculation required to directly invert the Hessian matrix is ​​very large (for example, if you want to accurately solve a two-dimensional model of size nx×nz, you need to perform nx×nz nx×nz The current computing hardware still cannot meet this demand. Therefore, an effective optimization scheme is needed to obtain low-cost and high-precision approximate Hessian operator. Operto proposed the Gauss-Newton method to approximate the Hessian operator, but this method has high requirements for computation and storage. Limited Memory (BFGS) uses gradients to construct approximate Hessian operators during iterations. However, this method consumes a lot of memory and is not conducive to processing large-scale data. In order to save memory consumption, L-BFGS estimates the approximate Hessian operator through the model and gradient information of the most recent few times. However, the matrix obtained by this method is positive definite, but the correct Hessian operator is non-positive definite. For highly ill-conditioned FWI, this positive definite approximation may lead to unsatisfactory inversion results. Unlike the above methods, the truncated Newton method (TRN) does not need to explicitly solve the Hessian operator. It only needs to iteratively update the inexact (truncated) solution related to the Newton equation through a small number of conjugate gradients in each iteration of FWI. The approximate Hessian operator can be constructed at a relatively low cost. Summary of the invention

[0005] In order to reconstruct the velocity and Q parameters with higher accuracy, the present invention proposes a viscoacoustic full waveform inversion method based on velocity-attenuation decoupling. A dual-parameter QFWI for simultaneous inversion of velocity and Q is proposed for the fractional-order viscoacoustic equation; the truncated Newton method is introduced to effectively approximate the Hessian operator, alleviate the problem of strong coupling effect of velocity and Q, and improve the inversion accuracy of dual parameters; the judgment condition is introduced into the TRN algorithm to avoid insufficient or excessive updates in the descending direction; this method can effectively suppress crosstalk, and can better restore the velocity and attenuation model, which can meet the requirements of the oil and gas industry for underground medium parameters.

[0006] The technical solution adopted by the present invention is: a viscoacoustic full waveform inversion method based on velocity-attenuation decoupling, the viscoacoustic full waveform inversion method comprising the following steps:

[0007] Step 1: Set the earthquake source, detection point, earthquake record information, and read the observation data;

[0008] Step 2: Calculate the fractional-order viscoacoustic equation in the forward direction according to time, obtain simulation data (Formula 2), and store part of the wave field; use Formula 9 to obtain the accompanying source, and then use Formula 7 to reversely calculate the accompanying wave field in time; use the partial wave field stored in Process 2 to reversely reconstruct the forward wave field and calculate Formula 5 and Formula 6; use Formula 4 to calculate the velocity and Q gradient; use Algorithm 1 to calculate the descent direction; use Formula 10 to update the model parameters;

[0009] Obtaining the objective function gradient, for QFWI, the least squares objective function can be expressed as:

[0010]

[0011] Where m represents the model parameter in space X; d and p represent the observed and simulated data wave fields respectively, and x r Represents the detector position;

[0012] The simulated wave field p can be obtained by forward modeling of seismic waves, and the viscoacoustic wave equation with fractional-order Laplace operator can be expressed as:

[0013] Ap=f (2)

[0014]

[0015] In the formula, f represents the earthquake source, c represents the medium velocity, ω0 is the reference angular frequency, and operator A is a partial differential equation (PDEs) describing the propagation of waves underground. The velocity and Q in this equation can be expressed explicitly, and it contains fewer parameters, which is conducive to the calculation of the gradient. Using the Lagrange multiplier method, the gradient G based on the objective function (1) can be expressed as:

[0016]

[0017] In the formula, <·,·> represents the inner product, μ represents the accompanying wave field, and for different parameters, the second term in the inner product can be expressed as:

[0018]

[0019]

[0020] The companion wave field of the first term in the inner product can be solved by the following formula:

[0021]

[0022] In the formula, Indicates accompanying, f * Indicates the accompanying earthquake source, represents the adjoint operator, which determines the scattering mode of the model parameters. Different wave equation systems and different characterization parameters have different scattering modes. The adjoint operator in this method can be expressed as:

[0023]

[0024] From Formula 7, we can see that the expression of the associated earthquake source is related to the form of the objective function, which can be expressed as:

[0025] f * =pd (9)

[0026] Since QFWI can be viewed as a PDE constrained optimization problem, the updated model parameter m can be expressed as:

[0027] m k+1 =m k +a k d k (10)

[0028] Where, d k and α k They represent the descent direction and descent step length at the kth iteration respectively. The step length can be estimated by linear search, and the model update amount d k It can be solved by the gradient method (Formula 11) or the Newton method (Formula 12):

[0029] d k =-P(m k )G(m k ) (11)

[0030] d k =-H -1 (m k )G(m k ) (12)

[0031] In the formula, P represents preconditioning, G represents gradient, and H represents Hessian matrix;

[0032] The gradient method (Formula 11) that ignores the effect of the Hessian matrix has a relatively low computational complexity, but the number of iterations required to reach convergence is very high. Although preconditioning methods such as L-BFGS can be used to accelerate the convergence speed, this preconditioning has limited effect on the decoupling of multi-parameter crosstalk. The Newton method (Formula 12) can alleviate problems such as data band limitation, uneven illumination, and multi-parameter coupling, and it can improve the convergence of the inversion. However, most of these methods require the inverse of the Hessian matrix, which requires a huge amount of computation, and current computing hardware still cannot meet this requirement.

[0033] Therefore, in order to solve the crosstalk problem of simultaneous inversion of multiple parameters in QFWI, the second key technology of this method is to use TRN to approximate the Newton descent direction. This method can avoid directly obtaining the inverse of the Hessian matrix without a matrix. The flowchart of using TRN to obtain the descent direction is shown in Algorithm 1.

[0034] The algorithm is:

[0035]

[0036] The algorithm gradually obtains the Newton descent direction in the internal conjugate gradient loop by introducing the product of the Hessian operator and the vector λ (Hessian vector product, Hλ); using the second-order adjoint state method, the Hessian vector product Hλ can be expressed as:

[0037]

[0038]

[0039]

[0040] In the formula, Ω, denote the computational domain and model parameters respectively, H(m)λ denotes the Hessian vector product of the m parameter. In order to ensure the convergence of the algorithm and avoid over-solving the Newton descent direction, a truncation condition is set:

[0041] ||H(m)λ+G(m)||>0.9||G(m)|| (16)

[0042] In addition, when encountering curvature associated with negative eigenvalues ​​of the operator, that is, when α2≤0, the inner conjugate gradient loop must stop, and the descent direction is the result of the last inner conjugate gradient. The above truncation condition ensures that the model update calculated by the inner conjugate gradient loop is always guaranteed to be in the descent direction of the objective function;

[0043] Step 3: If the maximum number of iterations is not reached, repeat step 2 until the maximum number of iterations is reached and the reconstructed speed and Q are obtained.

[0044] The beneficial effects of the present invention are as follows: a viscoacoustic full waveform inversion method based on velocity-attenuation decoupling is proposed. Different from the traditional acoustic wave FWI, the algorithm of the present method uses a constant fractional-order wave equation as the wave equation system, and can simultaneously reconstruct the velocity and attenuation parameters (Q); in order to avoid the crosstalk problem of simultaneous inversion of multiple parameters, a local optimization strategy based on TRN is designed, which can avoid directly solving the Hessian operator, but efficiently calculate the Hessian vector product in a matrix-free manner by using the second-order adjoint state method, and obtain the descent direction of the model parameters by using the Newton-conjugate gradient framework, and introduce an effective conjugate gradient termination condition to avoid over-update of the descent direction; the velocity and Q parameters with high resolution can be restored; the decoupling technology based on TRN can effectively alleviate the coupling effect in the simultaneous inversion of multiple parameters and improve the inversion accuracy; the truncation condition in the TRN algorithm can automatically determine whether the accurate descent direction is obtained, avoid over-update or insufficient update, and reduce the complexity of the algorithm; the present method can effectively suppress crosstalk, and can better restore the velocity and attenuation model, and can meet the requirements of the oil and gas industry for underground medium parameters. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1It is a flow chart of the viscoacoustic full waveform inversion method based on velocity-decay decoupling;

[0046] Figure 2 It is the velocity gradient map of the difference method for gradient verification;

[0047] Figure 3 It is the velocity gradient diagram of the adjoint method for gradient verification;

[0048] Figure 4 is the velocity gradient residual map for gradient verification;

[0049] Figure 5 is the gradient verified by the difference method Q gradient;

[0050] Figure 6 It is the Q gradient graph of the adjoint method for gradient verification;

[0051] Figure 7 is the Q gradient residual of gradient verification;

[0052] Figure 8 It is the velocity diagram of the difference method verified by the Hessian vector product;

[0053] Fig. 9 It is the adjoint method velocity diagram verified by the Hessian vector product;

[0054] Fig.10 It is the velocity residual map verified by the Hessian vector product;

[0055] Fig.11 It is the difference method Q diagram verified by Hessian vector product;

[0056] Fig.12 It is the adjoint method Q diagram verified by Hessian vector product;

[0057] Fig.13 is the Q residual graph verified by the Hessian vector product;

[0058] Fig.14 It is the real velocity model diagram;

[0059] Fig.15 It is the real Q model diagram;

[0060] Fig.16 is the initial velocity model diagram;

[0061] Fig.17 is the initial Q model diagram;

[0062] Fig.18 It is based on the velocity direction map of the gradient method;

[0063] Fig.19 It is based on the gradient method Q-direction map;

[0064] Fig. 20is the TRN decoupled velocity pattern;

[0065] Fig.21 is the TRN decoupled Q pattern;

[0066] Fig. 22 This is the velocity inversion result diagram based on TRN;

[0067] Fig.23 This is the Q inversion result diagram based on TRN;

[0068] Fig.24 This is the result map of velocity inversion based on gradient;

[0069] Fig.25 This is the gradient-based Q inversion result diagram. DETAILED DESCRIPTION

[0070] Embodiment 1

[0071] Reference Figure 1 , a viscoacoustic full waveform inversion method based on velocity-attenuation decoupling, the viscoacoustic full waveform inversion method comprising the following steps:

[0072] Step 1: Set the earthquake source, detection point, earthquake record information, and read the observation data;

[0073] Step 2: Calculate the fractional-order viscoacoustic equation in the forward direction according to time, obtain simulation data (Formula 2), and store part of the wave field; use Formula 9 to obtain the accompanying source, and then use Formula 7 to reversely calculate the accompanying wave field in time; use the partial wave field stored in Process 2 to reversely reconstruct the forward wave field and calculate Formula 5 and Formula 6; use Formula 4 to calculate the velocity and Q gradient; use Algorithm 1 to calculate the descent direction; use Formula 10 to update the model parameters;

[0074] Obtaining the objective function gradient, for QFWI, the least squares objective function can be expressed as:

[0075]

[0076] Where m represents the model parameter in space X; d and p represent the observed and simulated data wave fields respectively, and x r Represents the detector position;

[0077] The simulated wave field p can be obtained by forward modeling of seismic waves, and the viscoacoustic wave equation with fractional-order Laplace operator can be expressed as:

[0078] Ap=f (2)

[0079]

[0080] In the formula, f represents the earthquake source, c represents the medium velocity, ω0 is the reference angular frequency, and operator A is a partial differential equation (PDEs) describing the propagation of waves underground. The velocity and Q in this equation can be expressed explicitly, and it contains fewer parameters, which is conducive to the calculation of the gradient. Using the Lagrange multiplier method, the gradient G based on the objective function (1) can be expressed as:

[0081]

[0082] In the formula, <·,·> represents the inner product, μ represents the accompanying wave field, and for different parameters, the second term in the inner product can be expressed as:

[0083]

[0084]

[0085] The companion wave field of the first term in the inner product can be solved by the following formula:

[0086]

[0087] In the formula, Indicates accompanying, f * Indicates the accompanying earthquake source, represents the adjoint operator, which determines the scattering mode of the model parameters. Different wave equation systems and different characterization parameters have different scattering modes. The adjoint operator in this method can be expressed as:

[0088]

[0089] From Formula 7, we can see that the expression of the associated earthquake source is related to the form of the objective function, which can be expressed as:

[0090] f * =pd (9)

[0091] Since QFWI can be viewed as a PDE constrained optimization problem, the updated model parameter m can be expressed as:

[0092] m k+1 =m k +a k d k (10)

[0093] Where, d k and α k They represent the descent direction and descent step length at the kth iteration respectively. The step length can be estimated by linear search, and the model update amount d k It can be solved by the gradient method (Formula 11) or the Newton method (Formula 12):

[0094] d k =-P(mk )G(m k ) (11)

[0095] d k =-H -1 (m k )G(m k ) (12)

[0096] In the formula, P represents preconditioning, G represents gradient, and H represents Hessian matrix;

[0097] The gradient method (Formula 11) that ignores the effect of the Hessian matrix has a relatively low computational complexity, but the number of iterations required to reach convergence is very high. Although preconditioning methods such as L-BFGS can be used to accelerate the convergence speed, this preconditioning has limited effect on the decoupling of multi-parameter crosstalk. The Newton method (Formula 12) can alleviate problems such as data band limitation, uneven illumination, and multi-parameter coupling, and it can improve the convergence of the inversion. However, most of these methods require the inverse of the Hessian matrix, which requires a huge amount of computation, and current computing hardware still cannot meet this requirement.

[0098] Therefore, in order to solve the crosstalk problem of simultaneous inversion of multiple parameters in QFWI, the second key technology of this method is to use TRN to approximate the Newton descent direction. This method can avoid directly obtaining the inverse of the Hessian matrix without a matrix. The flowchart of using TRN to obtain the descent direction is shown in Algorithm 1.

[0099] The algorithm is:

[0100]

[0101] The algorithm gradually obtains the Newton descent direction in the internal conjugate gradient loop by introducing the product of the Hessian operator and the vector λ (Hessian vector product, Hλ); using the second-order adjoint state method, the Hessian vector product Hλ can be expressed as:

[0102]

[0103]

[0104]

[0105] In the formula, Ω, denote the computational domain and model parameters, respectively. H(m)λ denotes the Hessian vector product of the m parameter. In order to ensure the convergence of the algorithm and avoid over-solving the Newton descent direction, a truncation condition is set:

[0106] ||H(m)λ+G(m)||>0.9||G(m)|| (16)

[0107] In addition, when encountering curvature associated with negative eigenvalues ​​of the operator, that is, when α2≤0, the inner conjugate gradient loop must stop, and the descent direction is the result of the last inner conjugate gradient. The above truncation condition ensures that the model update calculated by the inner conjugate gradient loop is always guaranteed to be in the descent direction of the objective function;

[0108] Step 3: If the maximum number of iterations is not reached, repeat step 2 until the maximum number of iterations is reached and the reconstructed speed and Q are obtained.

[0109] Embodiment 2

[0110] In order to verify the correctness of the gradient and Hessian operator calculations in this method, a uniform velocity model is designed for verification, with the horizontal and vertical sizes and grid spacing of 300m and 10m respectively. A Ricker wavelet with a main frequency of 15Hz is used as the excitation sound source, the time sampling rate is 1ms, the true velocity and initial velocity are 2050m / s and 2000m / s, and the true and initial Q are 80 and 100 respectively. 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, formula 17):

[0111]

[0112] In the formula, It represents the gradient of the i-th point obtained by the definition method. Therefore, the definition method requires point-by-point wave field simulation to obtain the gradient of all points, while the adjoint state method only requires 2 wave field simulations to calculate the gradient of all points. i represents the disturbance parameter.

[0113] In this experiment, the disturbance velocity and disturbance Q are set to 1m / s and 20 respectively. In order to observe the first Fresnel zone between the source and the detector, the sound source and the detector are placed at (x=50m, y=50m) and (x=250m, y=250m) of the model respectively. Figure 2-7 As shown in the figure, both are displayed on the same scale. This shows that the gradients generated by the two methods are almost the same, and the velocity gradient is abnormal only at the source location, so the residual of the velocity gradient can be ignored. Since the perturbation of Q is relatively large, there is a small residual in the Q gradient. However, if the perturbation of Q is further reduced, the gradient energy of Q will be submerged in the numerical noise. However, as Figure 5 and Figure 6 As shown, the similarity is sufficient to verify the effectiveness of our adjoint state method for gradient calculation. From the definition of the Hessian matrix, we know that the Hessian matrix is ​​equal to the derivative of the gradient:

[0114]

[0115] In the formula, H i,j(m) represents the Hessian matrix element at position (i, j) in the model. When the vector λ in H(m)λ=(0,…,1,…0), H(m)λ≈H i，j (m).

[0116] The verification of Hessian vector product is as follows Figure 8-13 As shown in , the disturbance point is set at the center of the model. Figure 8-10 It can be seen that the result of the Hessian vector product of the adjoint state method is highly consistent with the difference method, with only two abnormal points at the source and the wave builder positions; while the Hessian vector product of Q has a disturbance anomaly at the disturbance point. But it is enough to prove the accuracy of the Hessian vector product.

[0117] Embodiment 3

[0118] In the simultaneous inversion of multiple parameters, the first problem to be solved is the crosstalk between different parameters. We verify the decoupling effect of Hessian vector product in multi-parameter inversion through anomaly model. The model size is 150×150 grid points, and the spatial spacing is 10m. The velocity model and Q model are as follows Figure 14-17 shown. Figure 14-17 Represent the real velocity model, the real Q model, the initial velocity model and the initial Q model respectively. A Ricker wavelet with a main frequency of 15 Hz is used as the excitation sound source, and the time sampling rate is 1 ms. 15 shots are evenly arranged on the surface, and the detectors are arranged around the model. The descent direction of the first iteration is as follows Figure 18-21 As shown, Fig.18 and Fig.19 It is a gradient-based method. Fig. 20 and Fig.21 It is based on the TRN method. Fig.18 and Fig.19 It can be seen that in the descent direction of the gradient method, there are different degrees of crosstalk between different parameters. Although the anomaly in the velocity descent direction is relatively accurate, in the attenuation descent direction, the velocity anomaly is incorrectly mapped to the attenuation direction, and even the energy far exceeds the attenuation anomaly. If this direction is used directly for parameter update, the Q parameter may not be correctly inverted. Fig. 20 and Fig.21 It can be observed that the energy of the anomalies with different parameters is more focused, the crosstalk problem between different parameters is effectively alleviated, and the attenuation anomalies are accurately located. The inversion results obtained in different descending directions (iteration 10 times) are shown in Figure 2. Figure 22-25 As shown, Figure 22-25 They represent the velocity inversion results based on TRN, the Q inversion results based on TRN, the velocity inversion results based on gradient, and the Q inversion results based on gradient. Figure 22-25 It can be seen that the velocity parameters reconstructed by the gradient-based method contain a weak attenuation disturbance ( Fig.24 ),and Fig.25 The attenuation parameters in the above method contain a large number of velocity "footprints", which seriously interfere with the correct reconstruction of the attenuation parameters. In this method, both velocity and attenuation anomalies are well reconstructed, and there is no obvious crosstalk problem ( Fig. 22 and Fig.23 ). Therefore, this method can effectively suppress the crosstalk between different parameters in QFWI and can better restore the velocity model and attenuation model.

Claims

1. A viscoacoustic full waveform inversion method based on velocity-attenuation decoupling, characterized by: The viscoacoustic full waveform inversion method comprises the following steps: Step 1: Set the earthquake source, detection point, earthquake record information, and read the observation data; Step 2: Solve the fractional-order viscoacoustic equation forward in time to obtain simulation data; obtain the accompanying wave field; construct the velocity and Q gradient; obtain the descent direction and update the model parameters; Obtaining the objective function gradient, for QFWI, the least squares objective function can be expressed as: Where m represents the model parameter in space X; d and p represent the observed and simulated data wave fields respectively, and x r Represents the detector position; The simulated wave field p can be obtained by forward modeling of seismic waves, and the viscoacoustic wave equation with fractional-order Laplace operator can be expressed as: Ap=f (2) In the formula, f represents the earthquake source, c represents the medium velocity, ω0 is the reference angular frequency, operator A is the partial differential equation describing the wave propagation in the underground. Using the Lagrange multiplier method, the gradient G based on the objective function (1) can be expressed as: In the formula, <·,·> represents the inner product, μ represents the accompanying wave field, and for different parameters, the second term in the inner product can be expressed as: The companion wave field of the first term in the inner product can be solved by the following formula: In the formula, Indicates accompanying, f * Indicates the accompanying earthquake source, represents the adjoint operator. In this method, the adjoint operator can be expressed as: The associated earthquake source is expressed as: f * =p-d (9) The updated model parameter m can be expressed as: m k+1 =m k +a k d k (10) Where, d k and α k They represent the descent direction and descent step length at the kth iteration respectively. The step length can be estimated by linear search, and the descent direction is obtained using TRN; The algorithm is: By introducing the product of the Hessian operator and the vector λ, the Newton descent direction is gradually obtained in the internal conjugate gradient loop; using the second-order adjoint state method, the Hessian vector product Hλ can be expressed as: In the formula, Ω, denote the computational domain and model parameters respectively, H(m)λ denotes the Hessian vector product of the m parameter. In order to ensure the convergence of the algorithm and avoid over-solving the Newton descent direction, a truncation condition is set: ||H(m)λ+G(m)||>0.9||G(m)|| (16); Step 3: If the maximum number of iterations is not reached, repeat step 2 until the maximum number of iterations is reached and the reconstructed speed and Q are obtained.

Images