An optimized and accelerated least squares reverse time migration imaging method

Abstract

The invention discloses an optimized accelerated least squares reverse time migration imaging method, which relates to the technical field of petroleum geophysical exploration, and includes a least squares reverse time migration imaging method based on a safe Anderson acceleration type I and a least squares reverse time migration imaging method based on an Anderson acceleration type II. The least squares reverse time migration imaging method based on a safe Anderson acceleration type I, combined with Powell regularization, restart check and safety assurance steps, greatly accelerates the convergence speed of data residuals; the least squares reverse time migration imaging method based on an Anderson acceleration type II, inputs an initial reflectivity model, memory parameters and fixed point operators, uses the initial reflectivity model and the fixed point operator to calculate the initial gradient, finds the minimum value of the number of iterations and the memory parameter in each iteration, saves the gradient vector of the minimum value number in each iteration, then corrects the gradient of each iteration, and uses the corrected gradient to update the reflectivity of each iteration for imaging.

Description

Technical Field

[0001] The invention relates to the technical field of petroleum geophysical exploration, and in particular to a least squares reverse time migration imaging method based on optimization acceleration. Background Art

[0002] With the development of exploration technology, geological targets and geological conditions have become more complex, and exploration data have become huge. Although least squares reverse time migration has become mature and has many applications in production, the calculation amount of least squares reverse time migration is large and it takes a lot of time to process data.

[0003] Traditional least squares reverse time migration methods usually use the steepest descent method or the conjugate gradient method to calculate the gradient. Although these methods can update the gradient and image the underground structure, they also have disadvantages, such as large amount of calculation. Each iteration needs to calculate the iteration step, which requires a lot of calculation cost and a long calculation time; the residual convergence speed is slow. The conjugate gradient method can only consider the effect of the previous iteration in each iteration, without considering the previous iterations, resulting in a slow convergence speed of the iterative residual.

[0004] Therefore, it is urgent to find a least squares reverse time migration imaging method that can solve the above problems. Summary of the invention

[0005] In order to solve the above technical problems, the present invention discloses an optimized accelerated least squares reverse time migration imaging method.

[0006] To achieve the above object, the present invention adopts the following technical solutions:

[0007] A least squares reverse time migration imaging method based on optimization acceleration comprises the following steps:

[0008] The acoustic wave equation for least squares reverse time migration in the frequency domain is:

[0009]

[0010] Among them, δ(xx η ) and δ(zz η ) are the impulse functions in the x and z directions respectively, is the slowness, ω is the angular frequency, W(x,z,ω) is the source function, and U(x,z,ω) is the wave field;

[0011] Equation (1) is discretized and the Cartesian coordinate system at 0° and 45° is used to differentiate the approximate Laplace term as follows:

[0012]

[0013] Among them, a is the weighting factor, is the Laplace operator;

[0014] The finite difference method of optimizing nine points is used to calculate the difference grid. The center point weight of the finite difference method of optimizing nine points is:

[0015] U i,j =cU i,j +d(U i,j-1 +U i,j+1 +U i-1,j +U i+1,j )+e(U i-1,j-1 +U i+1,j-1 +U i+1,j+1 +U i-1,j+1 )(3)

[0016] Among them, c, d, and e are all weighted factors, and the relationship between c, d, and e is

[0017] c+4d+4e=1 (4)

[0018] Substituting equations (2), (3) and (4) into equation (1), the discrete frequency domain least squares reverse time migration acoustic wave equation of the optimized nine-point finite difference method is obtained as follows:

[0019]

[0020] Among them, ρ, μ x , μ z , σ x , σ z All are optimized nine-point finite difference coefficients;

[0021] After adding the perfect matching layer, the optimized nine-point finite difference coefficients are expressed as:

[0022]

[0023] Among them, the values ​​of the optimized nine-point finite difference weighting factors are a=0.5461, c=0.6248, d=0.09381, ζ x and z is the attenuation function, △x and △z are the spacing of the computational grid in the x and z directions respectively;

[0024] The objective error function of the least squares reverse time migration method is defined as:

[0025]

[0026] Where L is the forward modeling operator, d obs (x,z,ω) is the observed data;

[0027] The model error is:

[0028]

[0029] The model error is calculated using a least squares reverse time migration imaging method based on improved Anderson acceleration.

[0030] Optionally, the least squares reverse time migration imaging method based on improved Anderson acceleration includes a least squares reverse time migration imaging method based on safe Anderson acceleration type I and a least squares reverse time migration imaging method based on Anderson acceleration type II.

[0031] Anderson acceleration type 1 is mainly achieved by reducing ||B k -I|| F To approximate g, it is subject to the following constraints:

[0032] B k *R k =Y k .

[0033] Optionally, the least squares reverse time migration imaging method based on the safe Anderson acceleration type 1 comprises the following steps:

[0034] Select the initial reflectivity model m0, memory parameter M≥1, fixed point operator f, set m1=f(m0),M k =min(M,k), k is the number of iterations, g χ =m χ -f(m χ ),y χ =g χ+1 -g χ ,r χ =m χ+1 -m χ , where χ = kM k ,…,k-1;

[0035] B k The solution is:

[0036]

[0037] The update strategy for the next iteration is

[0038] m k+1 =m k -(B k ) -1 *g k (10)

[0039] From equation (9), we can obtain (B k ) -1,as follows:

[0040]

[0041] Equation (10) can be rewritten as:

[0042]

[0043] in,

[0044] definition from Start Update as follows:

[0045]

[0046] in, Using Powell's method to define θ k ,as follows:

[0047]

[0048] in,

[0049] Restart the inspection strategy and set the counting parameter M A , set M A Increase by 1 each iteration until M = M A +1 or

[0050]

[0051] When the above conditions are met, reset M A To 0, set and H k =I;

[0052] Safety assurance step, calculate the second norm of g0 Then check whether the residual norm of each iteration is small enough; N A For counting parameters, if the following conditions are met

[0053]

[0054] N A Increase by 1, the next iteration becomes

[0055]

[0056] If condition (16) is not satisfied, the β equilibrium function is used to update m k+1 ,as follows

[0057]

[0058] Optionally, the anti-drift equation is obtained by an adjoint state method:

[0059]

[0060] Among them, U0(x,z,w) is the background wave field, U1(x,z,ω) is the first-order scattered wave field, m(x,z) is the reflectivity model, and d cal (x, ω) is the calculated data, and U1(x, 1, ω) represents the shot record at the detection point position (z = 1).

[0061] Optionally, the initial objective function gradient is also included, expressed as:

[0062]

[0063] Optionally, the least squares reverse time migration imaging method based on Anderson acceleration type II comprises the following steps:

[0064] Select reflectivity model m, memory parameter M, m is a number with n x ×n z A vector of grid points;

[0065] Optimizing weight vector Where k is the number of iterations, M k =min(M,k) is the minimum value of M and k, is a vector, is a numerical value. This optimization weight vector is obtained by solving the constrained optimization problem as follows:

[0066]

[0067] in, For the xth iteration, g χ =m χ -f(m χ ), f(m χ ) is m χ Fixed-point iteration function of ;

[0068] Update reflectivity model m,m k+1 Defined as

[0069]

[0070] in, is the first (kM k +x) times the reflectivity;

[0071] By reducing the parameters, equation (19) is transformed into the following unconstrained optimization problem:

[0072]

[0073] in, is a column vector, is the matrix of the kth iteration, and for χ=kM k ,…,k-1, there are y χ =g χ+1 -g χ ;

[0074] Assuming that the matrix Y is full rank, the unconstrained optimization problem is solved by the following equation

[0075]

[0076] Through equation (22), the optimization weight is obtained as

[0077]

[0078] By comparing the relationship between the constrained optimization problem (19) and equation (20), we can use the unconstrained optimization problem (21) to obtain m k+1 for

[0079]

[0080] Since f(m χ )=m χ -g χ and χ =g χ+1 -g χ ,definition For every χ=kM k ,…,k-1, all have r χ =m χ+1 -m χ ;

[0081] Anderson acceleration type II is expressed as

[0082]

[0083] Equation (25) is modified as

[0084] m k+1 =m k -H k *g k (26)

[0085] in,

[0086] The beneficial effects of the present invention are:

[0087] 1. The present invention proposes a least squares reverse time migration imaging method based on the safe Anderson acceleration type 1. The safe Anderson acceleration type 1 method combines Powell regularization, restart check and safety assurance steps. Powell regularization ensures the non-singularity of Anderson acceleration type 1, and the restart check strategy ensures strong linear dependence. Both of these contribute to the stability of Anderson acceleration type 1. The safety assurance step checks the decrease of data residuals and accelerates convergence. The data residual convergence speed of the safe Anderson acceleration type 1 method of the present invention is about twice as fast as the traditional LBFGS method, which can greatly accelerate the convergence speed of data residuals and reduce the calculation amount of least squares reverse time migration, thereby improving the calculation efficiency of least squares reverse time migration, solving the problems of large calculation amount of least squares reverse time migration and unstable data residual oscillation of Anderson acceleration type 1, explaining the calculation cost, and helping to improve exploration production efficiency in practical applications;

[0088] 2. The least squares reverse time migration imaging method based on Anderson acceleration type II of the present invention firstly needs to input the initial reflectivity model, memory parameters and fixed point operators, then uses the initial reflectivity model and the fixed point operator to find the initial gradient, finds the minimum value of the number of iterations and the memory parameters in each iteration, saves the gradient vector of the minimum value in each iteration, and then uses the Anderson acceleration type II acceleration method to correct the gradient of each iteration, and then can use the corrected gradient to update the reflectivity of each iteration, so as to perform imaging. Traditional least squares reverse time migration gradient update methods such as the steepest descent method and the conjugate gradient method need to calculate the step length, and calculating the step length takes a lot of computing time. The Anderson acceleration type II method uses a fixed step length and only modifies the gradient, which can reduce the time of calculating the step length, so it can save a lot of computing, which can help greatly improve the efficiency of the least squares reverse time migration calculation, and can also obtain accurate high-resolution and high signal-to-noise ratio imaging results. BRIEF DESCRIPTION OF THE DRAWINGS

[0089] Figure 1 are (a) the velocity model and (b) the true reflectivity model of the Marmousi model in Example 1 of the present invention;

[0090] Figure 2 The fastest descent method in Example 1 of the present invention, 5μ SD1 Anderson Accelerator Type 1, 5μ SD1 Safe Anderson acceleration type 1 data residual convergence curve, where μ SD1 represents the step size of the first iteration of the steepest descent method, 5μ SD1 Indicates the step size of the first iteration of the 5-fold steepest descent method;

[0091] Figure 3 is a data residual convergence curve of the safe Anderson acceleration type 1 with different step lengths in Example 1 of the present invention;

[0092] Figure 4 It is the data residual convergence curve of the conjugate gradient method, LBFGS, and safe Anderson acceleration type 1 in Example 1 of the present invention;

[0093] Figure 5 It is the model error convergence curve of the conjugate gradient method, LBFGS, and safe Anderson acceleration type 1 in Example 1 of the present invention;

[0094] Figure 6 1 is a time histogram of the conjugate gradient method, LBFGS, and safe Anderson acceleration type 1 at 10 iterations (white) and 20 iterations (black) in Example 1 of the present invention;

[0095] Figure 7 It is (a) 20-iteration conjugate gradient method in Example 1 of the present invention; (a) 20-iteration LBFGS; (c) 10-iteration safe Anderson acceleration type 1; (d) 20-iteration safe Anderson acceleration type 1 least squares reverse time migration image;

[0096] Figure 8 are (a) the velocity model and (b) the true reflectivity model of the Marmousi model in Example 2 of the present invention;

[0097] Fig. 9 are (a) the shot record of the true reflectivity in Example 2 of the present invention; (b) the shot record using the Anderson Acceleration Type II method after 20 iterations; (c) the difference between the Anderson Acceleration Type II and true reflectivity shot records;

[0098] Fig.10 It is a graph showing the variation of the residuals of the steepest descent method, the conjugate gradient method and the Anderson acceleration type II with the number of iterations in Example 2 of the present invention;

[0099] Fig.11 This is a graph showing the variation of the model error of the steepest descent method, the conjugate gradient method and the Anderson acceleration type II with the number of iterations in Example 2 of the present invention;

[0100] Fig.12 1 is a time histogram of the steepest descent method, the conjugate gradient method, and the Anderson acceleration type II in Example 2 of the present invention at 10 iterations (white) and 20 iterations (black);

[0101] Fig.13 It is the least squares reverse time migration image after 20 iterations of (a) the steepest descent method; (b) the conjugate gradient method; and (c) the Anderson acceleration type II in Example 2 of the present invention. DETAILED DESCRIPTION

[0102] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0103] Example 1

[0104] An optimized and accelerated least squares reverse time migration imaging method comprises the following steps:

[0105] The acoustic wave equation for least squares reverse time migration in the frequency domain is:

[0106]

[0107] Among them, δ(xx η ) and δ(zz η ) are the impulse functions in the x and z directions respectively, is the slowness, ω is the angular frequency, W(x,z,ω) is the source function, and U(x,z,ω) is the wave field;

[0108] Equation (1) is discretized and the Cartesian coordinate system at 0° and 45° is used to differentiate the approximate Laplace term as follows:

[0109]

[0110] Among them, a is the weighting factor, is the Laplace operator;

[0111] The finite difference method of optimizing nine points is used to calculate the difference grid. The center point weight of the finite difference method of optimizing nine points is:

[0112] U i,j =cU i,j +d(U i,j-1 +U i,j+1 +U i-1,j +U i+1,j )+e(U i-1,j-1 +U i+1,j-1 +U i+1,j+1 +U i-1,j+1 )(3)

[0113] Among them, c, d, and e are all weighted factors, and the relationship between c, d, and e is

[0114] c+4d+4e=1 (4)

[0115] Substituting equations (2), (3) and (4) into equation (1), the discrete frequency domain least squares reverse time migration acoustic wave equation of the optimized nine-point finite difference method is obtained as follows:

[0116]

[0117] Among them, ρ, μ x , μ z , σ x , σ z All are optimized nine-point finite difference coefficients;

[0118] After adding the perfect matching layer, the optimized nine-point finite difference coefficients are expressed as:

[0119]

[0120] Among them, the values ​​of the optimized nine-point finite difference weighting factors are a=0.5461, c=0.6248, d=0.09381, ζ x and z is the attenuation function, △x and △z are the spacing of the computational grid in the x and z directions respectively;

[0121] The objective error function of the least squares reverse time migration method is defined as:

[0122]

[0123] Where L is the forward modeling operator, d obs (x,z,ω) is the observed data;

[0124] The model error is:

[0125]

[0126] The least squares reverse time migration imaging method based on the safe Anderson acceleration type 1 includes the following steps:

[0127] Least squares reverse time migration using the secure Anderson accelerated type 1 algorithm

[0128]

[0129]

[0130] Application Cases

[0131] This example uses the Marmousi model to compare the performance and efficiency of the safe Anderson acceleration type 1 least squares reverse time migration and other methods. The size of the Marmousi model is 3680m×2000m, the two-dimensional grid spacing is 10m, the time step is 2ms, the sampling points are 1500, the total sampling time is 3s, the source is a Ricker wavelet with a main frequency of 30Hz, and the sampling frequency spacing is The least squares reverse time migration experiment of the Marmousi model has 18 shots, each shot point spacing is 200m, there are 368 receivers, the receiver point spacing is 10m, and the boundary of the perfect matching layer is 500m. The safe Anderson acceleration type I memory parameter M=5, and the other parameters are τ=0.001, β=0.1, D=1000, ε=0.001. The computer used in this embodiment is a 22-core Intel(R) Xeon(R) CPU E5-2699 v4 processor, and the least squares reverse time migration in the frequency domain is implemented on Matlab software, where: Figure 1 (a) is the velocity model of the Marmousi model, Figure 1 (b) is the true reflectance model of the Marmousi model.

[0132] Figure 2 Demonstrated the steepest descent method, 5μ SD1 Anderson Accelerator Type 1, 5μ SD1 The residual convergence curve of the safe Anderson accelerated type I least squares reverse time migration data, where μ SD1 represents the step size of the first iteration of the steepest descent method, μ SD1 The number in front represents the multiple. It can be seen that the Anderson acceleration type 1 has unstable data residual oscillation phenomenon, while the data residual curve of the safe Anderson acceleration type 1 can steadily decrease, which proves that the safe Anderson acceleration type 1 can eliminate the data residual oscillation phenomenon and make the Anderson acceleration type 1 stable.

[0133] Figure 3 The safe Anderson acceleration type 1 data residual convergence curves for different step sizes are shown, from which we can see that the step size is 10μ SD1 The data residual convergence curve can drop faster.

[0134] Figure 4 and Figure 5 The data residual and model error convergence curves of the conjugate gradient method, LBFGS and safe Anderson type 1 are shown respectively. Both figures show that the safe Anderson type 1 converges faster than the other two methods, and 10 iterations of the safe Anderson type 1 are approximately equivalent to the effect of 20 iterations of LBFGS. Figure 4It mainly focuses on the difference between the forward simulated gun records and the real gun records in each iteration of different methods. Compared with the other two methods, the data residual of the safe Anderson acceleration type 1 decreases faster, indicating that by updating the gradient, the safe Anderson acceleration type 1 can obtain the correct offset image faster. Figure 5 Mainly aiming at the difference between different methods and the real reflectivity model after different number of iterations, a normalization method is adopted. The closer the model error curve is to zero, the more accurate the image obtained by least squares reverse time migration is. The model error curve of the safe Anderson acceleration type 1 decreases much faster than that of other methods, which can indicate that the reflectivity model imaged by the method of the present invention is more accurate.

[0135] Figure 6 It shows the time consumed by 10 and 20 iterations of the conjugate gradient method, LBFGS and safe Anderson acceleration type I. It can be seen that the safe Anderson acceleration type I consumes less time than the conjugate gradient method and LBFGS method, and saves about 1 / 3 of the calculation time compared with the conjugate gradient method. This is mainly because the safe Anderson acceleration type I uses a fixed step size and does not need to additionally calculate the step size for each iteration. Therefore, it saves the time for calculating the step size and improves the computational efficiency of the least squares reverse time migration.

[0136] Figure 7 The images shown are the images of the least squares reverse time migration of the conjugate gradient method and LBFGS with 20 iterations, and the safe Anderson accelerated type 1 with 10 and 20 iterations. It can be seen that the image of the safe Anderson accelerated type 1 with 20 iterations has better imaging effect than the other two methods, with higher resolution and signal-to-noise ratio, and is closer to Figure 1 (b) The true reflectance model, and the image of the safe Anderson acceleration type I with 10 iterations and the image of the LBFGS method with 20 iterations are close.

[0137] Therefore, considering the large amount of calculation of least squares reverse time migration and the unstable problem of data residual oscillation of Anderson acceleration type I, the least squares reverse time migration imaging method based on the safe Anderson acceleration type I proposed in the present invention can eliminate the unstable problem of data residual oscillation of Anderson acceleration type I and improve the calculation efficiency of least squares reverse time migration, explain the calculation cost, and help to improve exploration production efficiency in practical applications.

[0138] Example 2

[0139] A least squares reverse time migration imaging method based on improved Anderson acceleration includes the following steps:

[0140] The acoustic wave equation for least squares reverse time migration in the frequency domain is:

[0141]

[0142] Among them, δ(xx η ) and δ(zz η ) are the impulse functions in the x and z directions respectively, is the slowness, ω is the angular frequency, W(x,z,ω) is the source function, and U(x,z,ω) is the wave field;

[0143] Equation (1) is discretized and the Cartesian coordinate system at 0° and 45° is used to differentiate the approximate Laplace term as follows:

[0144]

[0145] Among them, a is the weighting factor, is the Laplace operator;

[0146] The finite difference method of optimizing nine points is used to calculate the difference grid. The center point weight of the finite difference method of optimizing nine points is:

[0147] U i,j =cU i,j +d(U i,j-1 +U i,j+1 +U i-1,j +U i+1,j )+e(U i-1,j-1 +U i+1,j-1 +U i+1,j+1 +U i-1,j+1 )(3)

[0148] Among them, c, d, and e are all weighted factors, and the relationship between c, d, and e is

[0149] c+4d+4e=1 (4)

[0150] Substituting equations (2), (3) and (4) into equation (1), the discrete frequency domain least squares reverse time migration acoustic wave equation of the optimized nine-point finite difference method is obtained as follows:

[0151]

[0152] Among them, ρ, μ x , μ z , σ x , σ z All are optimized nine-point finite difference coefficients;

[0153] After adding the perfect matching layer, the optimized nine-point finite difference coefficients are expressed as:

[0154]

[0155] Among them, the values ​​of the optimized nine-point finite difference weighting factors are a=0.5461, c=0.6248, d=0.09381, ζ x and z is the attenuation function, △x and △z are the spacing of the computational grid in the x and z directions respectively;

[0156] The de-migration equation is obtained by the adjoint state method:

[0157]

[0158] Among them, U0(x,z,w) is the background wave field, U1(x,z,ω) is the first-order scattered wave field, m(x,z) is the reflectivity model, and d cal (x, ω) is the calculated data, U1(x, 1, ω) represents the shot record at the detection point position (z = 1);

[0159] The objective error function of the least squares reverse time migration method is defined as:

[0160]

[0161] Where L is the forward modeling operator, d obs (x,z,ω) is the observed data;

[0162] The model error is:

[0163]

[0164] The initial objective function gradient is expressed as:

[0165]

[0166] The least squares reverse time migration imaging method based on Anderson acceleration type II includes the following steps:

[0167] Anderson-accelerated type-2 algorithm for least squares reverse-time migration gradient update

[0168]

[0169] Application Cases

[0170] This embodiment uses the Marmousi model to compare the steepest descent method, the conjugate gradient method, and the Anderson acceleration type II. The model size is 368×200, the grid spacing in the horizontal and vertical directions is 10m, the time step is 2ms, the total sampling time is 3s, the source is a Ricker wavelet with a main frequency of 30Hz, there are 1500 sampling frequency points, and the sampling frequency spacing is 1 / 3Hz. Each shot has 368 receiving points, the detection points are 20m apart, the size of the perfect matching layer boundary condition is 500m, and the memory parameter M is selected as 5. This case is implemented on Matlab, and the computer used is a 22-core Intel(R) Xeon(R) CPU E5-2699 v4 processor, where, Figure 8 (a) is the velocity model of the Marmousi model, Figure 8 (b) is the true reflectance model of the Marmousi model.

[0171] Fig. 9 The difference between the forward modeled shot records made by least squares reverse time migration based on Anderson acceleration type II and the forward modeled shot records made by true reflectivity is demonstrated. It can be seen that after 20 iterations using Anderson acceleration type II, the reverse-migrated shot records can be close to the shot records made by true reflectivity, which proves that the Anderson acceleration type II method of the present invention is correct.

[0172] Fig.10 and Fig.11 The residuals and model errors of the steepest descent method, conjugate gradient method, and Anderson accelerated type II are shown respectively. Both figures show that Anderson accelerated type II converges faster than the other two methods. Fig.10 The residual convergence diagram is mainly for the difference between the forward simulation shot record and the real shot record of each iteration. It can be seen from it that the residual of the Anderson acceleration type II method of the present invention can be reduced quickly, which means that the correct migration image can be obtained quickly. Fig.11 It can be seen from the model error diagram that the curve of Anderson acceleration type II decreases much faster than that of the other two methods, which can explain that the reflectivity model imaged by the method of the present invention is closer to the real reflectivity model.

[0173] Fig.12 It shows the time taken for 10 and 20 iterations of the steepest descent method, conjugate gradient method and Anderson type II acceleration. From the 20 iterations, it can be seen that the Anderson type II acceleration method saves about 1 / 3 of the calculation time compared with the steepest descent method and the conjugate gradient method. This is mainly because the Anderson type II acceleration does not need to calculate the step size of each iteration, thus saving the step size calculation time, which can greatly improve the efficiency of least squares reverse time migration.

[0174] Fig.13The images shown are the least squares reverse time migration images of the steepest descent method, the conjugate gradient method, and the Anderson acceleration type II with 20 iterations. It can be seen that the Anderson acceleration type II can achieve very good imaging effects after 20 iterations, which is basically close to Figure 8 The real reflectivity model of (b) also has better imaging effect than the other two methods, with higher signal-to-noise ratio and resolution.

[0175] Therefore, considering the defect of large computational complexity of least squares reverse time migration, the present invention proposes a least squares reverse time migration imaging method based on Anderson acceleration type II, which can greatly improve the computational efficiency of least squares reverse time migration, save computational costs for seismic exploration, and help improve production and exploration efficiency.

[0176] Of course, the above description is not a limitation of the present invention, and the present invention is not limited to the above examples. Changes, modifications, additions or substitutions made by technicians in this technical field within the essential scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. An optimized and accelerated least squares reverse time migration imaging method, characterized in that: The steps include: The acoustic wave equation for least squares reverse time migration in the frequency domain is: Among them, δ(xx ε ) and δ(zz ε ) are the impulse functions in the x and z directions respectively, is the slowness, ω is the angular frequency, W(x,z,ω) is the source function, and U(x,z,ω) is the wave field; Equation (1) is discretized and the Cartesian coordinate system at 0° and 45° is used to differentiate the approximate Laplace term as follows: Among them, a is the weighting factor, is the Laplace operator; The finite difference method of optimizing nine points is used to calculate the difference grid. The center point weight of the finite difference method of optimizing nine points is: IN i,j =cU i,j +d(U i,j-1 +U i,j+1 +U i-1,j +U i+1,j )+e(U i-1,j-1 +U i+1,j-1 +U i+1,j+1 +U i-1,j+1 ) (3) Among them, c, d, and e are all weighted factors, and the relationship between c, d, and e is c+4d+4e=1 (4) Substituting equations (2), (3) and (4) into equation (1), the discrete frequency domain least squares reverse time migration acoustic wave equation of the optimized nine-point finite difference method is obtained as follows: Among them, ρ, μ x , μ z , x , z All are optimized nine-point finite difference coefficients; After adding the perfect matching layer, the optimized nine-point finite difference coefficients are expressed as: Among them, the values ​​of the optimized nine-point finite difference weighting factors are a = 0.5461, c = 0.6248, d = 0.09381, δ x and δ z is the attenuation function, Δx and Δz are the spacing of the computational grid in the x and z directions respectively; The objective error function of the least squares reverse time migration method is defined as: Where L is the forward modeling operator, d obs (x,w) is the observed data; The model error is: The model error is calculated using a least squares reverse time migration imaging method based on improved Anderson acceleration.

2. The optimized accelerated least squares reverse time migration imaging method according to claim 1, characterized in that: The least squares reverse time migration imaging method based on improved Anderson acceleration includes a least squares reverse time migration imaging method based on safe Anderson acceleration type I and a least squares reverse time migration imaging method based on Anderson acceleration type II.

3. The optimized accelerated least squares reverse time migration imaging method according to claim 2, characterized in that: The least squares reverse time migration imaging method based on the safe Anderson acceleration type 1 comprises the following steps: Select the initial reflectivity model m0, memory parameter M≥1, fixed point operator f, set m1=f(m0),M k =min(M,k), k is the number of iterations, g χ =m χ -f(m χ ),y χ =g χ+1 -g χ ,r χ =m χ+1 -m χ , where χ = kM k ,…,k-1; B k The solution is: The update strategy for the next iteration is m k+1 =m k -(B k ) -1 *g k (10) From equation (9), we can obtain (B k ) -1 ,as follows: Equation (10) can be rewritten as: in, definition from Start Update as follows: in, Using the Powell approach to define ζ k ,as follows: in, Restart the inspection strategy and set the counting parameter M A , set M A Increase by 1 each iteration until M = M A +1 or When the above conditions are met, reset M A To 0, set and H k =I; N A For counting parameters, if the following conditions are met N A Increase by 1, the next iteration becomes If condition (16) is not satisfied, the β equilibrium function is used to update m k+1 ,as follows 4. The optimized accelerated least squares reverse time migration imaging method according to claim 3, characterized in that: It also includes the de-migration equation obtained by the adjoint state method: Among them, U0(x,z,w) is the background wave field, U1(x,z,ω) is the first-order scattered wave field, m(x,z) is the reflectivity model, and d cal (x, ω) is the calculated data, and U1(x, 1, ω) represents the shot record at the detection point position (z = 1).

5. The optimized and accelerated least squares reverse time migration imaging method according to claim 4, characterized in that: It also includes the initial objective function gradient, expressed as:

6. The optimized accelerated least squares reverse time migration imaging method according to claim 5, characterized in that: The least squares reverse time migration imaging method based on Anderson acceleration type II comprises the following steps: Select the initial reflectivity model m0, memory parameters, fixed point operator f, and set Among them, μ0 is the initial step size; Optimizing weight vector Where k is the number of iterations, M k =min(M,k) is the minimum value of M and k, is a vector, is a numerical value, and the optimization weight vector is obtained by solving the constrained optimization problem as follows: in, For the xth iteration, g χ =m χ -f(m χ ), f(m χ ) is m χ Fixed-point iteration function of ; Update reflectivity model m,m k+1 Defined as in, is the first (kM k +x) times the reflectivity; By reducing the parameters, equation (19) is transformed into the following unconstrained optimization problem: in, is a column vector, is the matrix of the kth iteration, and for χ=kM k ,…,k-1, there are y χ =g χ+1 -g χ ; Assuming that the matrix Y is full rank, the unconstrained optimization problem is solved by the following equation Through equation (22), the optimization weight is obtained as By comparing the relationship between the constrained optimization problem (19) and equation (20), we can use the unconstrained optimization problem (21) to obtain m k+1 for Since f(m χ )=m χ -g χ and χ =g χ+1 -g χ ,definition For every χ=kM k ,…,k-1, all have r χ =m χ+1 -m χ ; Anderson acceleration type II is expressed as Equation (25) is modified as m k+1 =m k -H k *g k (26) in,

Images