Nonlinear least squares reverse time migration method based on linear Bern approximation constraints

Abstract

The present invention provides a nonlinear least squares reverse time migration method based on linear Berne approximation constraints, comprising: step 1, reconstructing a first-order scattering wave field through observation data and a first-order scattering wave equation under a given velocity disturbance; step 2, deriving an update gradient for calculating the velocity disturbance based on the first-order scattering wave field; step 3, discussing the convexity of the inversion problem of the nonlinear LSRTM based on the linear Berne approximation constraints; step 4, obtaining a high-precision and high-amplitude-preserving seismic profile based on the velocity disturbance update gradient calculated according to the first-order scattering wave field expression. The nonlinear least squares reverse time migration method based on linear Berne approximation constraints ensures the accuracy of the reconstruction of the first-order scattering wave field, thereby making the velocity disturbance update gradient calculated according to the first-order scattering wave field expression more accurate, so that we can obtain a high-precision and high-amplitude-preserving seismic profile.

Description

Technical Field

[0001] The invention relates to the technical field of petroleum geophysical exploration, and in particular to a nonlinear least squares reverse time migration method based on linear Bern approximation constraints. Background Art

[0002] Least squares reverse time migration can give high-resolution and high-amplitude imaging results by minimizing the difference between simulated data and observed data, and it is less likely to produce low-frequency noise and artifacts than ordinary reverse time migration methods. The LSRTM theory is based on the first-order scattering wave assumption under the Born approximation, and its corresponding scattering wave expression can only describe the wave field under weak scattering potential and small scatterers. However, the actual propagation process of seismic waves is a nonlinear process, involving multiple wave phenomena, and the underground medium is also very complex, and the conditions of weak scattering potential and small scatterers are generally difficult to meet. When the full wave field information is used for inversion, the high-order scattering waves generated by strong scattering potential and large scatterers, as well as direct waves, return waves and other waves will produce artifacts with the first-order scattering waves, affecting the imaging quality.

[0003] In the Chinese patent application with application number: CN202010756625.6, an optimized least squares reverse time migration imaging method is involved, which is implemented in accordance with the following steps: Step 1, using the wave equation to migrate and reverse migrate the seismic data; Step 2, using the conjugate gradient method to construct a least squares framework, and introducing an optimization factor to improve the convergence speed of the conjugate gradient method; Step 3, using the GPU's CUDA computing architecture to further improve the computational efficiency of the optimized LSRTM imaging.

[0004] In the Chinese patent application with application number: CN201811051319.1, a least squares migration imaging method and device based on the L-BFGS algorithm are involved. The method includes: performing reverse time migration imaging according to the acoustic wave velocity model of the target exploration area and the common shot point record of the acoustic wave velocity model to obtain an initial imaging result; calculating the record residual of the synthetic seismic record and the original observation record of the target exploration area according to the initial imaging result; calculating the gradient of the target functional with respect to the initial imaging result according to the record residual; obtaining a second imaging result through the L-BFGS algorithm according to the initial imaging result and the gradient of the target functional with respect to the initial imaging result; and judging whether the record residual and the number of calculations of the L-BFGS algorithm meet the preset termination condition.

[0005] In the Chinese patent application with application number: CN201610855634.4, a least squares migration imaging optimization method and system are involved. The method includes: obtaining the seismic data structure correlation coefficient based on the seismic signal energy information measurement within the time window, the seismic signal contrast measurement within the time window, and the structure similarity measurement within the time window; based on the seismic data structure correlation coefficient, establishing a least squares migration imaging model weighted by the seismic data structure correlation coefficient; and solving the least squares migration imaging model weighted by the seismic data structure correlation coefficient to obtain seismic reflection coefficient imaging.

[0006] The above existing technologies are all significantly different from the present invention and fail to solve the technical problem we want to solve. Therefore, we have invented a new nonlinear least squares reverse time migration method based on linear Bern approximation constraints. Summary of the invention

[0007] The purpose of the present invention is to provide a nonlinear least squares reverse time migration method based on linear Berne approximation constraints, which is used to obtain a high-precision and high-amplitude-preserving seismic profile.

[0008] The object of the present invention can be achieved by the following technical measures: a nonlinear least squares reverse time migration method based on a linear Berne approximation constraint, the nonlinear least squares reverse time migration method based on a linear Berne approximation constraint comprises:

[0009] Step 1: Under the given velocity disturbance, reconstruct the first-order scattered wave field through the observed data and the first-order scattered wave equation;

[0010] Step 2, deriving an updated gradient of velocity disturbance calculation based on the first-order scattered wave field;

[0011] Step 3, discuss the convexity of the inversion problem of nonlinear least squares reverse time migration (LSRTM) based on linear Bern approximation constraints;

[0012] Step 4: Update the gradient of the velocity perturbation calculated according to the first-order scattered wave field expression to obtain a high-precision and high-amplitude preserved seismic profile.

[0013] The purpose of the present invention can also be achieved by the following technical measures:

[0014] In step 1, the homogeneous Helmholtz equation corresponding to the propagation of the scalar field in the underground medium is:

[0015]

[0016] in, is the coordinate position, ω is the angular frequency, Δ is the Hamiltonian operator, u is the wave field value, Γ is the forward operator, and v is the velocity value.

[0017] Introducing background velocity field get:

[0018]

[0019] in, is the coordinate position, ω is the angular frequency, Δ is the Hamiltonian operator, u is the wave field value, v is the velocity value, and v 0 is the background speed.

[0020] The scattering potential is defined as:

[0021]

[0022] in, is the coordinate position, ω is the angular frequency, v is the velocity value, V is the scattering potential, v 0 is the background speed, This is the velocity disturbance of interest to LSRTM.

[0023] According to the scattering potential, the above formula is simplified to:

[0024]

[0025] in, is the coordinate position, ω is the angular frequency, Δ is the Hamiltonian operator, u is the wave field value, v 0 is the background velocity field, and V is the scattering potential.

[0026] It is characterized in that, in step 1, it is assumed that the incident wave field propagating in the background medium satisfies the homogeneous Helmholtz equation:

[0027]

[0028] in, is the coordinate position, ω is the angular frequency, Δ is the Hamiltonian operator, u is the wave field value, v 0 is the background velocity field.

[0029] Then equation (5) is simplified to a non-homogeneous Helmholtz equation, which describes the scattered wave field generated after the incident wave field propagating in the background medium encounters the scatterer:

[0030]

[0031] in, is the coordinate position, ω is the angular frequency, Δ is the Hamiltonian operator, u is the wave field value, v 0 is the background velocity field, u sis the scattered wave field, and V is the scattering potential.

[0032] Define Green's function:

[0033]

[0034] in, is the coordinate position, ω is the angular frequency, Δ is the Hamiltonian operator, G is the Green function, δ is the pulse function, is the moving coordinate position.

[0035] The analytical solution of the scattered wave field is obtained using Green's function:

[0036]

[0037] Among them, u s is the scattered wave field, V is the scattering potential, u is the wave field value, G is the Green function, is the coordinate position, is the moving coordinate position, ω is the angular frequency.

[0038] Since the above formula is an expression of a volume integral, the scattered wave field is controlled by the scattering potential and the size of the scatterer.

[0039] In step 1, under the assumption of the first-order Born approximation, that is, the scattering potential is weak and the scatterer is small, the generated scattered wave field is much smaller than the background wave field: The scattered wave field under the first-order Born approximation is obtained:

[0040]

[0041] in, is the first-order scattered wave field, V is the scattering potential, u is the wave field value, G is the Green function, is the moving coordinate position, ω is the angular frequency, is the coordinate position of the detection point, and Ω is the integral volume.

[0042] The corresponding full wave field is:

[0043]

[0044] Among them, u is the wave field value, u 0 is the background wave field value, V is the scattering potential, G is the Green function, ω is the angular frequency, is the coordinate position of the detection point, Ω is the integral volume, is the moving coordinate position.

[0045] Substituting formula (10) into formula (8), we get the second-order scattered wave expression. Repeating the above process can get the multi-order scattered wave field expression, and simply express the N-order scattered wave field in summation form:

[0046]

[0047] Among them, u s is the scattered wave field, is the i-th order scattered wave field, N is the total order, is the coordinate position, ω is the angular frequency.

[0048] In step 2, when the scattered wave equation under the first-order Born approximation is used as the forward operator, the LSRTM objective functional under the L-2 norm is expressed as:

[0049]

[0050] Where J is the target functional, is the first-order simulated scattered wave field, is the first-order scattered wave field, d n is the nth-order background wave field.

[0051] The corresponding velocity perturbation updated gradient is expressed as:

[0052]

[0053] in, is the velocity perturbation, J is the target functional, is the coordinate position, ω is the angular frequency, is the first-order simulated scattered wave field, is the first-order simulated wave field, d n is the nth-order background wave field, u 0 is the background wave field value, V is the scattering potential, G is the Green's function, is the coordinate position of the detection point.

[0054] Analysis shows that the above gradient consists of two parts: the gradient formed by the residual between the simulated background wave field and the simulated first-order scattered wave field and the first-order scattered wave field in the observed data:

[0055]

[0056] in, is the first-order gradient, is the coordinate position, ω is the angular frequency, u 0 is the background wave field value, G is the Green function, is the coordinate position of the detection point, is the moving coordinate position, is the first-order simulated scattered wave field, It is a first-order simulated wave field.

[0057] The gradient formed by the simulated background wave field, the background wave field in the observed data, and the observed data under the high-order Born approximation is:

[0058]

[0059] Among them, g n is the nth order gradient, ω is the angular frequency, u 0 is the background wave field value, is the coordinate position, is the coordinate position of the detection point, d n is the nth-order background wave field, G is the Green function, is the moving coordinate position.

[0060] The first term of the gradient describes the velocity perturbation corresponding to the first-order scattering wave described by the forward operator, while the second term appears as imaging noise under the first-order scattering wave assumption, that is, under the first-order Born approximation, other components of the observed data cannot be accurately described, and thus interfere with the velocity perturbation gradient together with the background wave field.

[0061] In step 2, in order to suppress the noise generated by the background wave field and the high-order scattered wave field in the observed data, the linear Bern approximation equation is added as a constraint term to the objective function of the nonlinear LSRTM to strengthen the constraint of the first-order scattered wave in the nonlinear LSRTM to obtain a high-precision, high-resolution profile.

[0062] In step 2, the improved LSRTM objective function is expressed as:

[0063]

[0064] Among them, λ is the constraint factor, which is adjusted according to the intensity of the background wave field and high-order scattered wave field in the observed data, J is the target functional, is the coordinate position, ω is the angular frequency, is the first-order simulated scattered wave field, is the first-order simulated wave field, d n is the nth-order background wave field, u 0 is the background wave field value, V is the scattering potential, G is the Green's function, is the coordinate position of the detection point, is the moving coordinate position.

[0065] After adding constraints, LSRTM is decomposed into two inversion problems: first, the first-order scattering wave field is inverted to obtain a relatively accurate first-order scattering wave field; then, on this basis, the updated gradient of the velocity perturbation is calculated; let the derivative of formula (16) with respect to the first-order scattering wave field be 0, and the reconstructed expression of the first-order scattering wave field can be obtained:

[0066]

[0067] Among them, λ is the constraint factor, which is adjusted according to the intensity of the background wave field and high-order scattered wave field in the observed data, J is the target functional, is the coordinate position, ω is the angular frequency, is the first-order simulated scattered wave field, is the first-order simulated wave field, d n is the nth-order background wave field, u 0 is the background wave field value, V is the scattering potential, G is the Green's function, is the coordinate position of the detection point, is the moving coordinate position.

[0068]

[0069] Among them, λ is the constraint factor, is the coordinate position, ω is the angular frequency, u s 1 is the first-order simulated scattered wave field, is the first-order simulated wave field, d n is the nth-order background wave field, u 0 is the background wave field value, V is the scattering potential, G is the Green's function, is the moving coordinate position.

[0070] From the analysis of formula (18), it can be seen that the reconstructed first-order scattered wave field is mainly composed of two parts: the wave field composed of observed data and the simulated first-order scattered wave field. The proportion of the two is controlled by the constraint factor: when the constraint factor is larger, the influence of the first-order scattered wave on the reconstructed scattered wave field becomes greater, thereby weakening the influence of the background wave field and the high-order scattered wave on the first-order reconstructed scattered wave.

[0071] In step 2, after accurately reconstructing the scattered wave field, the gradient of the objective function with respect to the velocity disturbance is obtained:

[0072]

[0073] Among them, λ is the constraint factor, J is the target functional, is the coordinate position, ω is the angular frequency, u 0 is the background wave field value, V is the scattering potential, is the moving coordinate position, G is the Green's function, is the velocity disturbance, is the first-order simulated scattering wave field, u is the wave field value, u 0 is the background wave field value.

[0074] Then the velocity disturbance is updated:

[0075]

[0076] Where α is the update step size, J is the target functional, is the speed disturbance.

[0077] In step 3, for the general inversion problem, the property of the second-order derivative of the objective function with respect to the inversion parameters, namely the Hessian matrix, determines the convexity of the inversion problem: when the Hessian matrix is ​​semi-positive definite, namely the objective function is fully convex, and the global extreme value of the objective function is found by the local optimization method; otherwise, namely the objective function is locally convex, the inversion problem is affected by the local extreme value and it is difficult to converge to the global extreme value.

[0078] In step 3, the Hessian matrix of the first-order scattered wave field based on the nonlinear LSRTM constrained by the linear Bern approximation is:

[0079]

[0080] Among them, λ is the constraint factor, J is the target functional, is the coordinate position, is the coordinate position, ω is the angular frequency, u 0 is the background wave field value, V is the scattering potential, G is the Green's function, is the moving coordinate position, is the first-order simulated scattered wave field, is the first-order simulated wave field, d n is the nth-order background wave field, is the coordinate position of the detection point.

[0081] From the expression of the first-order scattered wave field, it can be seen that the wave field is affected by the model parameters, scattering potential and observation system; the wave fields under different observation systems are unrelated, then:

[0082]

[0083] Among them, H is the second-order gradient, r is the detection point position, λ is the constraint factor, J is the target functional, is the coordinate position, is the coordinate position, ω is the angular frequency, is the first-order simulated scattered wave field.

[0084] That is, the nonlinear LSRTM based on the linear Born approximation constraint shows strong convexity for the first-order scattering wave field, and a relatively accurate first-order scattering wave field is obtained by using the local optimization method; in addition, the first-order scattering wave field under the Born approximation is linear with the velocity perturbation, so an accurate velocity perturbation update gradient is obtained based on the reconstructed first-order scattering wave field.

[0085] The nonlinear least squares reverse time migration method based on linear Berne approximation constraints in the present invention proposes to add the linear Berne approximation as a constraint term to the objective function of the nonlinear LSRTM, and suppress the illusion caused by the weak scattering potential and the small scatterer assumption by strengthening the constraints of the linear Berne approximation during the inversion process. The inversion method based on the linear Berne approximation constraints decomposes the nonlinear LSRTM problem into two inversion problems: first, the optimal solution is found in the first-order scattering wave field space, that is, under the condition of given velocity perturbation, an accurate first-order scattering wave field is reconstructed as much as possible through the observed data and the first-order scattering wave equation. Then, the updated gradient of the velocity perturbation is calculated based on the first-order scattering wave field. The new objective function shows strong convexity to the Hessian matrix of the first-order scattering wave field, which ensures the accuracy of the reconstruction of the first-order scattering wave field, so that the velocity perturbation update gradient calculated according to the first-order scattering wave field expression is also more accurate, so that we can obtain a high-precision and high-amplitude seismic profile.

[0086] The beneficial effects of the present invention are as follows: the present invention provides a nonlinear least squares reverse time migration method based on linear Bern approximation constraints, and reconstructs an accurate first-order scattering wave field as much as possible through observation data and the first-order scattering wave equation. Then, the updated gradient of the velocity disturbance is calculated based on the first-order scattering wave field. The new objective function shows strong convexity to the Hessian matrix of the first-order scattering wave field, which ensures the accuracy of the reconstruction of the first-order scattering wave field, thereby making the updated gradient of the velocity disturbance calculated according to the first-order scattering wave field expression more accurate, so that we can obtain a high-precision and high-amplitude-preserving seismic profile. BRIEF DESCRIPTION OF THE DRAWINGS

[0087] Figure 1 A schematic diagram of the Sigsbee2A model in a specific embodiment of the present invention;

[0088] Figure 2 It is a schematic diagram of the shot recording for inversion in the specific embodiment 1 of the present invention;

[0089] Figure 3 Schematic diagram of the Sigsbee2A migration velocity model in the specific embodiment 1 of the present invention;

[0090] Figure 4 It is a schematic diagram of the LSRTM result after 10 iterations in the specific embodiment 1 of the present invention;

[0091] Figure 5 For the specific embodiment 1 of the present invention Figure 4 A magnified view of the part above the salt;

[0092] Figure 6 For the specific embodiment 1 of the present invention Figure 4 A magnified view of the subsalt section;

[0093] Figure 7 It is a schematic diagram of the shot recording for inversion in the specific embodiment 2 of the present invention;

[0094] Figure 8 It is a schematic diagram of the Sigsbee2A migration velocity model in the specific embodiment 2 of the present invention;

[0095] Fig. 9 It is a schematic diagram of the LSRTM result after 10 iterations in the specific embodiment 2 of the present invention;

[0096] Fig.10 In the specific embodiment 2 of the present invention, Figure 4 A magnified view of the part above the salt;

[0097] Fig.11 In the specific embodiment 2 of the present invention, Figure 4 A magnified view of the subsalt section;

[0098] Fig.12 It is a schematic diagram of the shot recording for inversion in the specific embodiment 3 of the present invention;

[0099] Fig.13 It is a schematic diagram of the Sigsbee2A migration velocity model in the specific embodiment 3 of the present invention;

[0100] Fig.14 It is a schematic diagram of the LSRTM result after 10 iterations in the specific embodiment 3 of the present invention;

[0101] Fig.15 For the specific embodiment 3 of the present invention Figure 4 A magnified view of the part above the salt;

[0102] Fig.16 For the specific embodiment 3 of the present invention Figure 4 A magnified view of the subsalt section;

[0103] Fig.17 The flowchart is a specific embodiment of the nonlinear least squares reverse time migration method based on linear Berne approximation constraints of the present invention. DETAILED DESCRIPTION

[0104] It should be noted that the following detailed descriptions are exemplary and are intended to provide further explanation of the present invention. Unless otherwise specified, all technical and scientific terms used herein have the same meanings as those commonly understood by those skilled in the art to which the present invention belongs.

[0105] It should be noted that the terms used herein are only for describing specific embodiments and are not intended to limit exemplary embodiments according to the present invention. As used herein, unless the context clearly indicates otherwise, the singular form is also intended to include the plural form. In addition, it should be understood that when the terms "comprising" and / or "including" are used in this specification, it indicates the presence of features, steps, operations and / or combinations thereof.

[0106] like Fig.17 As shown, Fig.17 The flowchart of the nonlinear least squares reverse time migration method based on the linear Berne approximation constraint of the present invention is as follows:

[0107] Step 101, under a given velocity disturbance, reconstructing an accurate first-order scattered wave field as far as possible through observation data and a first-order scattered wave equation;

[0108] The homogeneous Helmholtz equation corresponding to the propagation of the scalar field in the underground medium is:

[0109]

[0110] in, is the coordinate position, ω is the angular frequency, Δ is the Hamiltonian operator, u is the wave field value, Γ is the forward operator, and v is the velocity value.

[0111] Introducing background velocity field We can get:

[0112]

[0113] in, is the coordinate position, ω is the angular frequency, Δ is the Hamiltonian operator, u is the wave field value, v is the velocity value, and v 0 is the background speed.

[0114] The scattering potential is defined as:

[0115]

[0116] in, is the coordinate position, ω is the angular frequency, v is the velocity value, V is the scattering potential, v 0 is the background speed, This is the velocity perturbation that LSRTM is interested in. According to the scattering potential, the above formula can be simplified as:

[0117]

[0118] in, is the coordinate position, ω is the angular frequency, Δ is the Hamiltonian operator, u is the wave field value, v 0 is the background velocity field, and V is the scattering potential.

[0119] Assume that the incident wave field propagating in the background medium satisfies the homogeneous Helmholtz equation:

[0120]

[0121] in, is the coordinate position, ω is the angular frequency, Δ is the Hamiltonian operator, u is the wave field value, v 0 is the background velocity field.

[0122] Then equation (5) is simplified to a non-homogeneous Helmholtz equation, which describes the scattered wave field generated after the incident wave field propagating in the background medium encounters the scatterer:

[0123]

[0124] in, is the coordinate position, ω is the angular frequency, Δ is the Hamiltonian operator, u is the wave field value, v 0 is the background velocity field, u s is the scattered wave field, and V is the scattering potential.

[0125] Define Green's function:

[0126]

[0127] in, is the coordinate position, ω is the angular frequency, Δ is the Hamiltonian operator, G is the Green function, δ is the pulse function, is the moving coordinate position.

[0128] The analytical solution of the scattered wave field can be obtained using Green's function:

[0129]

[0130] Among them, u s is the scattered wave field, V is the scattering potential, u is the wave field value, G is the Green function, is the coordinate position, is the moving coordinate position, and ω is the angular frequency.

[0131] Since the above formula is a volume integral expression, the scattered wave field is controlled by the scattering potential and the size of the scatterer. Under the assumption of the first-order Born approximation, that is, the scattering potential is weak and the scatterer is small, the generated scattered wave field is much smaller than the background wave field: We can obtain the scattered wave field under the first-order Born approximation:

[0132]

[0133] in, is the first-order scattered wave field, V is the scattering potential, u is the wave field value, G is the Green function, is the moving coordinate position, ω is the angular frequency, is the coordinate position of the detection point, and Ω is the integral volume.

[0134] The corresponding full wave field is:

[0135]

[0136] Among them, u is the wave field value, u 0 is the background wave field value, V is the scattering potential, G is the Green function, ω is the angular frequency, is the coordinate position of the detection point, Ω is the integral volume, is the moving coordinate position.

[0137] Substituting formula (10) into formula (8), we can get the expression of the second-order scattered wave. Repeating the above process can get the expression of multi-order scattered wave field, and simply express the N-order scattered wave field in summation form:

[0138]

[0139] Among them, u s is the scattered wave field, is the i-th order scattered wave field, N is the total order, is the coordinate position, ω is the angular frequency.

[0140] Step 102, deriving an update gradient of velocity disturbance calculation based on the first-order scattered wave field;

[0141] The updated gradient of the velocity disturbance is derived based on the first-order scattered wave field, including:

[0142] When the scattered wave equation under the first-order Born approximation is used as the forward operator, the LSRTM objective functional under the L-2 norm can be expressed as:

[0143]

[0144] Where J is the target functional, is the first-order simulated scattered wave field, is the first-order scattered wave field, d n is the nth-order background wave field.

[0145] The corresponding velocity perturbation updated gradient can be expressed as:

[0146]

[0147] in, is the velocity perturbation, J is the target functional, is the coordinate position, ω is the angular frequency, is the first-order simulated scattered wave field, is the first-order simulated wave field, d n is the nth-order background wave field, u 0 is the background wave field value, V is the scattering potential, G is the Green's function, is the coordinate position of the detection point.

[0148] Analysis shows that the above gradient consists of two parts: the gradient formed by the residual between the simulated background wave field and the simulated first-order scattered wave field and the first-order scattered wave field in the observed data:

[0149]

[0150] in, is the first-order gradient, is the coordinate position, ω is the angular frequency, u 0 is the background wave field value, G is the Green function, is the coordinate position of the detection point, is the moving coordinate position, is the first-order simulated scattered wave field, It is a first-order simulated wave field.

[0151] The gradient formed by the simulated background wave field, the background wave field in the observed data, and the observed data under the high-order Born approximation is:

[0152]

[0153] Among them, g n is the nth order gradient, ω is the angular frequency, u 0 is the background wave field value, is the coordinate position, is the coordinate position of the detection point, d n is the nth-order background wave field, G is the Green function, is the moving coordinate position.

[0154] Obviously, the first term of the gradient describes the velocity perturbation corresponding to the first-order scattered wave described by the forward operator, while the second term appears as imaging noise under the assumption of the first-order scattered wave, that is, under the first-order Born approximation, other components of the observed data cannot be accurately described, and the background wave field interferes with the velocity perturbation gradient. In order to suppress the noise generated by the background wave field and the high-order scattered wave field in the observed data, the linear Born approximation equation is added as a constraint term to the objective function of the nonlinear LSRTM, and the constraint of the first-order scattered wave in the nonlinear LSRTM is strengthened to obtain a high-precision, high-resolution profile. The improved LSRTM objective function can be expressed as:

[0155]

[0156] Among them, λ is the constraint factor, which is adjusted according to the intensity of the background wave field and high-order scattered wave field in the observed data, J is the target functional, is the coordinate position, ω is the angular frequency, is the first-order simulated scattered wave field, is the first-order simulated wave field, d n is the nth-order background wave field, u 0 is the background wave field value, V is the scattering potential, G is the Green's function, is the coordinate position of the detection point, is the moving coordinate position.

[0157] After adding constraints, LSRTM can be decomposed into two inversion problems: first, the first-order scattering wave field is inverted to obtain a relatively accurate first-order scattering wave field; then, the updated gradient of the velocity perturbation is calculated on this basis. Let the derivative of formula (16) for the first-order scattering wave field be 0, and the reconstruction expression of the first-order scattering wave field can be obtained:

[0158]

[0159] Among them, λ is the constraint factor, which is adjusted according to the intensity of the background wave field and high-order scattered wave field in the observed data, J is the target functional, is the coordinate position, ω is the angular frequency, is the first-order simulated scattered wave field, is the first-order simulated wave field, d n is the nth-order background wave field, u 0 is the background wave field value, V is the scattering potential, G is the Green's function, is the coordinate position of the detection point, is the moving coordinate position.

[0160]

[0161] Among them, λ is the constraint factor, is the coordinate position, ω is the angular frequency, is the first-order simulated scattered wave field, is the first-order simulated wave field, d n is the nth-order background wave field, u 0 is the background wave field value, V is the scattering potential, G is the Green's function, is the moving coordinate position.

[0162] From the analysis of formula (18), it can be seen that the reconstructed first-order scattering wave field is mainly composed of two parts: the wave field composed of the observed data and the simulated first-order scattering wave field. The proportion of the two is controlled by the constraint factor: when the constraint factor is larger, the influence of the first-order scattering wave on the reconstructed scattering wave field becomes greater, thereby weakening the influence of the background wave field and the high-order scattering wave on the first-order reconstructed scattering wave.

[0163] After accurately reconstructing the scattered wave field, the gradient of the objective function to the velocity disturbance can be obtained:

[0164]

[0165] Among them, λ is the constraint factor, J is the target functional, is the coordinate position, ω is the angular frequency, u 0 is the background wave field value, V is the scattering potential, is the moving coordinate position, G is the Green's function, is the velocity disturbance, is the first-order simulated scattering wave field, u is the wave field value, u 0 is the background wave field value.

[0166] Then the velocity disturbance can be updated:

[0167]

[0168] Where α is the update step size, J is the target functional, is the speed disturbance.

[0169] Step 103, discuss the convexity of the inversion problem of the nonlinear LSRTM based on the linear Bern approximation constraint; for general inversion problems, the properties of the second-order derivative (Hessian matrix) of the objective function with respect to the inversion parameters generally determine the convexity of the inversion problem: when the Hessian matrix is ​​positive semi-definite, that is, the objective function is fully convex, and the global extremum of the objective function can be found by the local optimization method; on the contrary, that is, the objective function is locally convex, the inversion problem is affected by the local extremum and it is difficult to converge to the global extremum. The Hessian matrix of the nonlinear LSRTM based on the linear Bern approximation constraint for the first-order scattering wave field is:

[0170]

[0171] Among them, λ is the constraint factor, J is the target functional, is the coordinate position, is the coordinate position, ω is the angular frequency, u 0 is the background wave field value, V is the scattering potential, G is the Green's function, is the moving coordinate position, is the first-order simulated scattered wave field, is the first-order simulated wave field, d n is the nth-order background wave field, is the coordinate position of the detection point.

[0172] From the expression of the first-order scattered wave field, it can be seen that the wave field is affected by the model parameters, scattering potential and observation system. The wave fields under different observation systems are unrelated, so:

[0173]

[0174] Among them, λ is the constraint factor, J is the target functional, is the coordinate position, is the coordinate position, ω is the angular frequency, is the first-order simulated scattered wave field.

[0175] That is, the nonlinear LSRTM based on the linear Born approximation constraint shows strong convexity for the first-order scattering wave field, and a relatively accurate first-order scattering wave field can be obtained by local optimization. In addition, the first-order scattering wave field under the Born approximation is linear with the velocity perturbation, so we can obtain an accurate velocity perturbation update gradient based on the reconstructed first-order scattering wave field.

[0176] Step 104, updating the gradient of the velocity disturbance calculated according to the first-order scattered wave field expression, and obtaining a high-precision and high-amplitude-preserving seismic profile.

[0177] Embodiment 1:

[0178] In the specific embodiment 1 of the present invention, the least squares reverse time migration based on the primary scattered wave constraint is first applied to the Sigsbee2A model ( Figure 1 ) to test its adaptability. The Sigsbee2A model contains a high-speed salt dome, so it cannot meet the assumptions of weak scattering potential and small scatterers, which can well illustrate the effectiveness of the method of the present invention. The model size is 401*1601, and the grid spacing is 21m vertically and 22m horizontally. The time domain acoustic wave equation is also used to synthesize seismic records, with a sampling time of 7s and a sampling interval of 1ms. 80 guns are placed on the surface, with a gun point interval of 440m, and full-aperture reception on the surface. Figure 2A shot record used for inversion is shown. This synthetic record contains multiple wave phenomena such as direct waves and return waves. Due to the existence of salt domes and the complexity of the deep model, the shot record also includes interlayer multiple waves, high-order scattered waves and other phenomena. The background velocity used for inversion is as follows Figure 3 As shown in the figure, the model has no obvious scatterers except for large salt domes. In order to reflect the importance of constraint factors to inversion, the inversion was performed using ordinary LSRTM, LSRTM with a constraint factor of 1000, and LSRTM with a constraint factor of 100000. The final inversion results are shown in Figure 1. Figure 4 As shown, Figure 4 (a) is the conventional LSRTM result; (b) is the LSRTM result based on the first-order scattered wave constraint (constraint factor is 1000); (c) is the LSRTM result based on the first-order scattered wave constraint (constraint factor is 100000). Also because the direct wave and the return wave energy are too strong, the inversion results of the deep layer are not obvious.

[0179] In order to more clearly compare the imaging results of the two methods. First, a part above the salt is intercepted for comparison: Due to the influence of direct waves and return waves, there is a lot of interference noise in the inversion results given by ordinary LSRTM, and the shallow arcing phenomenon is serious. After introducing the constraint of the first-order scattered wave, the interference noise is significantly reduced, and the resolution of the inversion result is further improved. With the further strengthening of the constraint of the first-order scattered wave, the noise generated by the background wave field such as the direct wave and the return wave has been basically eliminated, and the shallow layer position of the inversion result is clearer. Figure 5 For Figure 4 Magnified view of the part above the salt.

[0180] Then, some methods in the subsalt area were selected for precise comparison: due to the complexity of the observation data, the ordinary LSRTM results contain a lot of noise, and no scattering points are inverted. The inversion results of the LSRTM with added constraints have better fidelity and improved resolution, but no scattering points are inverted either. By further strengthening the constraints of the first-order scattering waves, LSRTM finally gave a better inversion result: the subsalt layer and scattering point positions are accurate, and the resolution is high. Figure 6 For Figure 4 Magnified view of the subsalt section.

[0181] Least squares reverse time migration can give high-resolution and high-amplitude imaging results by minimizing the difference between simulated data and observed data, and it is less likely to produce low-frequency noise and artifacts than ordinary reverse time migration methods. The LSRTM theory is based on the first-order scattering wave assumption under the Born approximation, and its corresponding scattering wave expression can only describe the wave field under weak scattering potential and small scatterers. However, the actual propagation process of seismic waves is a nonlinear process, involving multiple wave phenomena, and the underground medium is also very complex, and the conditions of weak scattering potential and small scatterers are generally difficult to meet. When the full wave field information is used for inversion, the high-order scattering waves generated by strong scattering potential and large scatterers, as well as direct waves, return waves and other waves will produce artifacts with the first-order scattering waves, affecting the imaging quality.

[0182] Therefore, taking into account the defects of the above-mentioned technologies, the present invention proposes to add the expression of the linear Berne approximation as a constraint term into the objective function of the nonlinear LSRTM, and to suppress the illusions caused by the weak scattering potential and the small scatterer assumption by enhancing the constraints of the first-order scattered waves during the inversion process.

[0183] Embodiment 2:

[0184] In the specific embodiment 2 of the present invention, the least squares reverse time migration based on the primary scattered wave constraint is first applied to the Sigsbee2A model ( Figure 1 ) to test its adaptability. The Sigsbee2A model contains a high-speed salt dome, so it cannot meet the assumptions of weak scattering potential and small scatterers, which can well illustrate the effectiveness of the method of the present invention. The model size is 201*801, and the grid spacing is 42m vertically and 44m horizontally. The time domain acoustic wave equation is also used to synthesize seismic records, with a sampling time of 7s and a sampling interval of 1ms. 80 guns are placed on the surface, with a gun point interval of 440m, and full-aperture reception on the surface. Figure 7 A shot record used for inversion is shown. This synthetic record contains multiple wave phenomena such as direct waves and return waves. Due to the existence of salt domes and the complexity of the deep model, the shot record also includes interlayer multiple waves, high-order scattered waves and other phenomena. The background velocity used for inversion is as follows Figure 8 As shown in the figure, the model has no obvious scatterers except for large salt domes. In order to reflect the importance of constraint factors to inversion, the inversion was performed using ordinary LSRTM, LSRTM with a constraint factor of 1000, and LSRTM with a constraint factor of 100000. The final inversion results are shown in Figure 1. Fig. 9 As shown, Fig. 9 (a) is the conventional LSRTM result; (b) is the LSRTM result based on the first-order scattered wave constraint (constraint factor is 1000); (c) is the LSRTM result based on the first-order scattered wave constraint (constraint factor is 100000). Also because the direct wave and the return wave energy are too strong, the inversion results of the deep layer are not obvious.

[0185] In order to more clearly compare the imaging results of the two methods. First, a part above the salt is intercepted for comparison: Due to the influence of direct waves and return waves, there is a lot of interference noise in the inversion results given by ordinary LSRTM, and the shallow arcing phenomenon is serious. After introducing the constraint of the first-order scattered wave, the interference noise is significantly reduced, and the resolution of the inversion result is further improved. With the further strengthening of the constraint of the first-order scattered wave, the noise generated by the background wave field such as the direct wave and the return wave has been basically eliminated, and the shallow layer position of the inversion result is clearer. Fig.12 For Fig. 9 Magnified view of the part above the salt.

[0186] Then, some methods in the subsalt area were selected for precise comparison: due to the complexity of the observation data, the ordinary LSRTM results contain a lot of noise, and no scattering points are inverted. The inversion results of the LSRTM with added constraints have better fidelity and improved resolution, but no scattering points are inverted either. By further strengthening the constraints of the first-order scattering waves, LSRTM finally gave a better inversion result: the subsalt layer and scattering point positions are accurate, and the resolution is high. Fig.11 For Fig. 9 Magnified view of the subsalt section.

[0187] Embodiment 3:

[0188] In the specific embodiment 3 of the present invention, the least squares reverse time migration based on the primary scattered wave constraint is first applied to the Sigsbee2A model ( Figure 1 ) to test its adaptability. The Sigsbee2A model contains a high-speed salt dome, so it cannot meet the assumptions of weak scattering potential and small scatterers, which can well illustrate the effectiveness of the method of the present invention. The model size is 101*401, and the grid spacing is 84m vertically and 88m horizontally. The time domain acoustic wave equation is also used to synthesize seismic records, with a sampling time of 7s and a sampling interval of 1ms. 80 guns are placed on the surface, with a gun point interval of 440m, and full-aperture reception on the surface. Figure 8 A shot record used for inversion is shown. This synthetic record contains multiple wave phenomena such as direct waves and return waves. Due to the existence of salt domes and the complexity of the deep model, the shot record also includes interlayer multiple waves, high-order scattered waves and other phenomena. The background velocity used for inversion is as follows Fig.14 As shown in the figure, the model has no obvious scatterers except for large salt domes. In order to reflect the importance of constraint factors to inversion, the inversion was performed using ordinary LSRTM, LSRTM with a constraint factor of 1000, and LSRTM with a constraint factor of 100000. The final inversion results are shown in Figure 1. Fig.14 As shown, Fig.14(a) is the conventional LSRTM result; (b) is the LSRTM result based on the first-order scattered wave constraint (constraint factor is 1000); (c) is the LSRTM result based on the first-order scattered wave constraint (constraint factor is 100000). Also because the direct wave and the return wave energy are too strong, the inversion results of the deep layer are not obvious.

[0189] In order to more clearly compare the imaging results of the two methods. First, a part above the salt is intercepted for comparison: Due to the influence of direct waves and return waves, there is a lot of interference noise in the inversion results given by ordinary LSRTM, and the shallow arcing phenomenon is serious. After introducing the constraint of the first-order scattered wave, the interference noise is significantly reduced, and the resolution of the inversion result is further improved. With the further strengthening of the constraint of the first-order scattered wave, the noise generated by the background wave field such as the direct wave and the return wave has been basically eliminated, and the shallow layer position of the inversion result is clearer. Fig.15 For Fig.14 Magnified view of the part above the salt.

[0190] Then, some methods in the subsalt area were selected for precise comparison: due to the complexity of the observation data, the ordinary LSRTM results contain a lot of noise, and no scattering points are inverted. The inversion results of the LSRTM with added constraints have better fidelity and improved resolution, but no scattering points are inverted either. By further strengthening the constraints of the first-order scattering waves, LSRTM finally gave a better inversion result: the subsalt layer and scattering point positions are accurate, and the resolution is high. Fig.16 For Fig.14 Magnified view of the subsalt section.

[0191] Finally, it should be noted that the above is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art can still modify the technical solutions recorded in the aforementioned embodiments or replace some of the technical features therein with equivalents. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

[0192] Except for the technical features described in the specification, all other technical features are known technologies to those skilled in the art.

Claims

1. A nonlinear least squares reverse time migration method based on linear Bern approximation constraints, characterized in that: The nonlinear least squares reverse time migration method based on linear Bern approximation constraints includes: Step 1: Under the given velocity disturbance, reconstruct the first-order scattered wave field through the observed data and the first-order scattered wave equation; Step 2, deriving an updated gradient of velocity disturbance calculation based on the first-order scattered wave field; Step 3, discuss the convexity of the inversion problem of nonlinear LSRTM based on the linear Bern approximation constraint; Step 4, updating the gradient of the velocity perturbation calculated according to the first-order scattered wave field expression to obtain a high-precision and high-amplitude-preserving seismic profile; In step 2, when the scattered wave equation under the first-order Born approximation is used as the forward operator, the LSRTM objective functional under the L-2 norm is expressed as: Where J is the target functional, is the first-order simulated scattered wave field, is the first-order scattered wave field, d n is the nth-order background wave field; The corresponding velocity perturbation updated gradient is expressed as: in, is the velocity perturbation, J is the target functional, is the coordinate position, ω is the angular frequency, is the first-order simulated scattered wave field, is the first-order simulated wave field, d n is the nth-order background wave field, u0 is the background wave field value, V is the scattering potential, G is the Green function, is the coordinate position of the detection point; Analysis shows that the above gradient consists of two parts: the gradient formed by the residual between the simulated background wave field and the simulated first-order scattered wave field and the first-order scattered wave field in the observed data: in, is the first-order gradient, is the coordinate position, ω is the angular frequency, u0 is the background wave, G is the Green function, is the coordinate position of the detection point, is the moving coordinate position, is the first-order simulated scattered wave field, is the first-order simulated wave field; The gradient formed by the simulated background wave field, the background wave field in the observed data, and the observed data under the high-order Born approximation is: Among them, g n is the nth order gradient, ω is the angular frequency, u0 is the background wave field value, is the coordinate position, is the coordinate position of the detection point, d n is the nth-order background wave field, G is the Green function, is the moving coordinate position; The first term of the gradient describes the velocity perturbation corresponding to the first-order scattered wave described by the forward operator, while the second term appears as imaging noise under the assumption of the first-order scattered wave, that is, under the first-order Born approximation, other components of the observed data cannot be accurately described, and thus interfere with the velocity perturbation gradient with the background wave field; In step 2, in order to suppress the noise generated by the background wave field and the high-order scattered wave field in the observed data, the linear Bern approximation equation is added as a constraint term to the objective function of the nonlinear LSRTM, and the constraint of the first-order scattered wave in the nonlinear LSRTM is strengthened to obtain a high-precision and high-resolution profile; In step 2, the improved LSRTM objective function is expressed as: Among them, λ is the constraint factor, which is adjusted according to the intensity of the background wave field and high-order scattered wave field in the observed data, J is the target functional, is the coordinate position, ω is the angular frequency, is the first-order simulated scattered wave field, is the first-order simulated wave field, d n is the nth order background wave field, u0 is the background wave field value, V is the scattering potential, G is the Green function, is the coordinate position of the detection point, is the moving coordinate position; After adding constraints, LSRTM is decomposed into two inversion problems: first, the first-order scattering wave field is inverted to obtain a relatively accurate first-order scattering wave field; then, on this basis, the updated gradient of the velocity perturbation is calculated; let the derivative of formula (16) with respect to the first-order scattering wave field be 0, and the reconstructed expression of the first-order scattering wave field can be obtained: Among them, λ is the constraint factor, which is adjusted according to the intensity of the background wave field and high-order scattered wave field in the observed data, J is the target functional, is the coordinate position, ω is the angular frequency, is the first-order simulated scattered wave field, is the first-order simulated wave field, d n is the nth order background wave field, u0 is the background wave field value, V is the scattering potential, G is the Green function, is the coordinate position of the detection point, is the moving coordinate position; Among them, λ is the constraint factor, is the coordinate position, ω is the angular frequency, is the first-order simulated scattered wave field, is the first-order simulated wave field, d n is the nth order background wave field, u0 is the background wave field value, V is the scattering potential, G is the Green function, is the moving coordinate position; From the analysis of formula (18), it can be seen that the reconstructed first-order scattered wave field is mainly composed of two parts: the wave field composed of observed data and the simulated first-order scattered wave field. The proportion of the two is controlled by the constraint factor: when the constraint factor is larger, the influence of the first-order scattered wave on the reconstructed scattered wave field becomes greater, thereby weakening the influence of the background wave field and the high-order scattered wave on the first-order reconstructed scattered wave.

2. The nonlinear least squares reverse time migration method based on linear Berne approximation constraints according to claim 1, characterized in that: In step 1, the homogeneous Helmholtz equation corresponding to the propagation of the scalar field in the underground medium is: in, is the coordinate position, ω is the angular frequency, Δ is the Hamiltonian operator, u is the wave field value, Γ is the forward operator, and v is the velocity value; Introducing background velocity field get: in, is the coordinate position, ω is the angular frequency, Δ is the Hamiltonian operator, u is the wave field value, v is the velocity value, and v0 is the background velocity; The scattering potential is defined as: in, is the coordinate position, ω is the angular frequency, v is the velocity value, V is the scattering potential, v0 is the background velocity, That is the velocity disturbance of interest to LSRTM; According to the scattering potential, the above formula is simplified to: in, is the coordinate position, ω is the angular frequency, Δ is the Hamiltonian operator, u is the wave field value, v0 is the background velocity field, and V is the scattering potential.

3. The nonlinear least squares reverse time migration method based on linear Berne approximation constraints according to claim 2, characterized in that: In step 1, it is assumed that the incident wave field propagating in the background medium satisfies the homogeneous Helmholtz equation: in, is the coordinate position, ω is the angular frequency, Δ is the Hamiltonian operator, u is the wave field value, and v0 is the background velocity field; Then equation (5) is simplified to a non-homogeneous Helmholtz equation, which describes the scattered wave field generated after the incident wave field propagating in the background medium encounters the scatterer: in, is the coordinate position, ω is the angular frequency, Δ is the Hamiltonian operator, u is the wave field value, v0 is the background velocity field, u s is the scattered wave field, V is the scattering potential; Define Green's function: in, is the coordinate position, ω is the angular frequency, Δ is the Hamiltonian operator, G is the Green function, δ is the pulse function, is the moving coordinate position; The analytical solution of the scattered wave field is obtained using Green's function: Among them, u s is the scattered wave field, V is the scattering potential, u is the wave field value, G is the Green function, is the coordinate position, is the moving coordinate position, ω is the angular frequency; Since the above formula is an expression of a volume integral, the scattered wave field is controlled by the scattering potential and the size of the scatterer.

4. The nonlinear least squares reverse time migration method based on linear Berne approximation constraints according to claim 3, characterized in that: In step 1, under the assumption of the first-order Born approximation, that is, the scattering potential is weak and the scatterer is small, the generated scattered wave field is much smaller than the background wave field: The scattered wave field under the first-order Born approximation is obtained: in, is the first-order scattered wave field, V is the scattering potential, u is the wave field value, G is the Green function, is the moving coordinate position, ω is the angular frequency, is the coordinate position of the detection point, Ω is the integral volume; The corresponding full wave field is: Among them, u is the wave field value, u0 is the background wave field value, V is the scattering potential, G is the Green function, ω is the angular frequency, is the coordinate position of the detection point, Ω is the integral volume, is the moving coordinate position; Substituting formula (10) into formula (8), we get the second-order scattered wave expression. Repeating the above process can get the multi-order scattered wave field expression, and simply express the N-order scattered wave field in summation form: Among them, u s is the scattered wave field, is the i-th order scattered wave field, N is the total order, is the coordinate position, ω is the angular frequency.

5. The nonlinear least squares reverse time migration method based on linear Berne approximation constraints according to claim 1, characterized in that: In step 2, after accurately reconstructing the scattered wave field, the gradient of the objective function with respect to the velocity disturbance is obtained: Among them, λ is the constraint factor, J is the target functional, is the coordinate position, ω is the angular frequency, V is the scattering potential, is the moving coordinate position, G is the Green function, is the velocity disturbance, is the first-order simulated scattering wave field, u is the wave field value, and u0 is the background wave field value; Then the velocity disturbance is updated: Among them, α is the update step size, J is the target functional, is the speed disturbance.

6. The nonlinear least squares reverse time migration method based on linear Berne approximation constraints according to claim 1, characterized in that: In step 3, for the inversion problem, the property of the second-order derivative of the objective function with respect to the inversion parameters, namely the Hessian matrix, determines the convexity of the inversion problem: when the Hessian matrix is ​​semi-positive definite, that is, the objective function is fully convex, and the global extreme value of the objective function is found by the local optimization method; otherwise, that is, the objective function is locally convex, the inversion problem is affected by the local extreme value and it is difficult to converge to the global extreme value.

7. The nonlinear least squares reverse time migration method based on linear Berne approximation constraints according to claim 6, characterized in that: In step 3, the Hessian matrix of the first-order scattered wave field based on the nonlinear LSRTM constrained by the linear Bern approximation is: Among them, λ is the constraint factor, J is the target functional, is the coordinate position, is the coordinate position, ω is the angular frequency, u0 is the background wave field value, V is the scattering potential, G is the Green's function, is the moving coordinate position, is the first-order simulated scattered wave field, is the first-order simulated wave field, d n is the nth-order background wave field, is the coordinate position of the detection point; From the expression of the first-order scattered wave field, it can be seen that the wave field is affected by the model parameters, scattering potential and observation system; the wave fields under different observation systems are unrelated, then: Among them, λ is the constraint factor, J is the target functional, is the coordinate position, is the coordinate position, ω is the angular frequency, is the first-order simulated scattering wave field, N grid is the grid point; That is, the nonlinear LSRTM based on the linear Born approximation constraint shows strong convexity for the first-order scattering wave field, and a relatively accurate first-order scattering wave field is obtained by using the local optimization method; in addition, the first-order scattering wave field under the Born approximation is linear with the velocity perturbation, so an accurate velocity perturbation update gradient is obtained based on the reconstructed first-order scattering wave field.

Images