CN115480294A - Nonlinear least square reverse time migration method based on linear Beren approximate constraint - Google Patents
Nonlinear least square reverse time migration method based on linear Beren approximate constraint Download PDFInfo
- Publication number
- CN115480294A CN115480294A CN202110602166.0A CN202110602166A CN115480294A CN 115480294 A CN115480294 A CN 115480294A CN 202110602166 A CN202110602166 A CN 202110602166A CN 115480294 A CN115480294 A CN 115480294A
- Authority
- CN
- China
- Prior art keywords
- order
- wave field
- scattering
- coordinate position
- background
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 56
- 230000005012 migration Effects 0.000 title claims abstract description 37
- 238000013508 migration Methods 0.000 title claims abstract description 37
- 238000004321 preservation Methods 0.000 claims abstract description 7
- 238000004613 tight binding model Methods 0.000 claims description 32
- 238000003384 imaging method Methods 0.000 claims description 24
- 238000004088 simulation Methods 0.000 claims description 16
- 239000011159 matrix material Substances 0.000 claims description 14
- 238000001514 detection method Methods 0.000 claims description 12
- 238000005457 optimization Methods 0.000 claims description 8
- 230000001902 propagating effect Effects 0.000 claims description 6
- 239000000523 sample Substances 0.000 claims description 4
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 3
- 230000000644 propagated effect Effects 0.000 claims 1
- 150000003839 salts Chemical class 0.000 description 30
- 238000010586 diagram Methods 0.000 description 7
- 238000005070 sampling Methods 0.000 description 6
- 238000004422 calculation algorithm Methods 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 3
- 230000002194 synthesizing effect Effects 0.000 description 3
- 239000002131 composite material Substances 0.000 description 2
- 238000002939 conjugate gradient method Methods 0.000 description 2
- 230000002708 enhancing effect Effects 0.000 description 2
- HPTJABJPZMULFH-UHFFFAOYSA-N 12-[(Cyclohexylcarbamoyl)amino]dodecanoic acid Chemical compound OC(=O)CCCCCCCCCCCNC(=O)NC1CCCCC1 HPTJABJPZMULFH-UHFFFAOYSA-N 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 238000011524 similarity measure Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/301—Analysis for determining seismic cross-sections or geostructures
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/303—Analysis for determining velocity profiles or travel times
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention provides a nonlinear least square reverse time migration method based on linear Bern approximate constraint, which comprises the following steps: step 1, reconstructing a first-order scattered wave field through observation data and a first-order scattered wave equation under the condition of given speed disturbance; step 2, deducing an updating gradient of the calculated speed disturbance based on the first-order scattering wave field; step 3, discussing the convexity of the inverse problem of the nonlinear LSRTM based on the linear Bern approximate constraint; and 4, updating the gradient according to the velocity disturbance calculated by the first-order scattering wave field expression to obtain a seismic section with high precision and high amplitude preservation. The nonlinear least square reverse time migration method based on the linear Berne approximate constraint ensures the accuracy of the reconstruction of the first-order scattering wave field, so that the velocity disturbance updating gradient calculated according to the first-order scattering wave field expression is more accurate, and a seismic section with high precision and high amplitude preservation can be obtained.
Description
Technical Field
The invention relates to the technical field of petroleum geophysical exploration, in particular to a nonlinear least square reverse time migration method based on linear Bern approximate constraint.
Background
Least squares reverse time migration can give high resolution and high amplitude preservation of imaging results by minimizing the difference between the simulated data and the observed data, and is less prone to low frequency noise and artifacts than conventional reverse time migration methods. The LSRTM theory is based on the assumption of first-order scattered waves under Born approximation, and its corresponding scattered wave expression can only describe weak scattering potentials and wave fields under small scatterers. However, the actual seismic wave propagation process is a nonlinear process, which includes various wave phenomena, and the underground medium is very complex, and the conditions of weak scattering potential and small scattering bodies are generally difficult to satisfy. When the full wave field information is used for inversion, high-order scattered waves generated by strong scattering potential and large scatterers, direct waves, rotating waves and other waves generate false images with first-order scattered waves, and imaging quality is affected.
In the application No.: chinese patent application CN202010756625.6 relates to an optimized least squares reverse time migration imaging method, which is specifically implemented according to the following steps: step 1, migration and reverse migration are carried out on seismic data by adopting a wave equation; step 2, constructing a least square frame by adopting a conjugate gradient method, and introducing an optimization factor to improve the convergence speed of the conjugate gradient method; and 3, further improving the calculation efficiency of the optimized LSRTM imaging by adopting a CUDA calculation framework of the GPU.
In the application No.: chinese patent application CN201811051319.1, which relates to a least square offset imaging method and device based on L-BFGS algorithm, the method includes: performing reverse time migration imaging according to the sound wave velocity model of the target exploration area and the common shot point record of the sound wave velocity model to obtain an initial imaging result; calculating the recording residual error 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 relative to the initial imaging result according to the recorded residual error; obtaining a second imaging result through an L-BFGS algorithm according to the initial imaging result and the gradient of the target functional about the initial imaging result; and judging whether the recorded residual error and the calculation times of the L-BFGS algorithm meet the preset termination condition.
In application No.: CN201610855634.4, chinese patent application, relates to a least squares offset imaging optimization method and system. The method comprises the following steps: obtaining a seismic data structure correlation coefficient based on the seismic signal energy information measure in the time window, the seismic signal contrast measure in the time window and the structure similarity measure in the time window; establishing a least squares migration imaging model weighted by the seismic data structure correlation coefficients based on the seismic data structure correlation coefficients; and solving a least square migration imaging model weighted by the seismic data structure correlation coefficient to obtain seismic reflection coefficient imaging.
The prior art is greatly different from the method, and the technical problem which is required to be solved by the inventor is not solved, so that the inventor invents a novel nonlinear least square reverse time migration method based on linear Bern approximate constraint.
Disclosure of Invention
The invention aims to provide a nonlinear least square reverse time migration method based on linear Bern approximate constraint, which is used for obtaining a high-precision and high-amplitude seismic section and is based on linear Bern approximate constraint.
The object of the invention can be achieved by the following technical measures: the nonlinear least square reverse time migration method based on the linear Bern approximate constraint comprises the following steps:
step 1, reconstructing a first-order scattered wave field through observation data and a first-order scattered wave equation under the condition of given speed disturbance;
step 2, deducing an updating gradient of the calculated speed disturbance based on the first-order scattering wave field;
and 4, updating the gradient according to the velocity disturbance calculated by the first-order scattering wave field expression to obtain a seismic section with high precision and high amplitude preservation.
The object of the invention can also be achieved by the following technical measures:
in step 1, the homogeneous Helmholtz equation for the scalar field of propagation in the subsurface medium is:
wherein,
is the coordinate position, omega is the angular frequency, delta is the hamiltonian, u is the wavefield value, Γ is the positive operator, v is the velocity value.
Introduction of background velocity field
Obtaining:
wherein,
is the coordinate position, omega is the angular frequency, delta is the Hamiltonian, u is the wavefield value, v is the velocity value 0 Is the background velocity.
The scattering potential is defined as:
wherein,
is the coordinate position, omega is the angular frequency, V is the velocity value, V is the scattering potential, V 0 In order to be the background velocity,
i.e. the velocity disturbance of interest to the LSRTM.
According to the scattering potential, the above formula is simplified as follows:
wherein,
is the coordinate position, omega is the angular frequency, delta is the Hamiltonian, u is the wave field value, v 0 For background velocity field, V is the scattering potential.
It is characterized in that in step 1, the incident wavefield propagating in the background medium is assumed to satisfy the homogeneous Helmholtz equation:
wherein,
is the coordinate position, ω is the angular frequency, Δ is the Hamiltonian, u is the wave field value, v 0 Is the background velocity field.
Equation (5) is reduced to a non-homogeneous Helmholtz equation describing the scattered wavefield that results after the incident wavefield propagating in the background medium encounters the scatterers:
wherein,
is the coordinate position, omega is the angular frequency, delta is the Hamiltonian, u is the wave field value, v 0 As background velocity field, u s Is the scattering wavefield and V is the scattering potential.
Define the Green function of Green:
wherein,
is the coordinate position, omega is the angular frequency, delta is the Hamiltonian, G is the Green function, delta is the impulse function,
is a mobile coordinate location.
Obtaining an analytic solution of the scattering wave field by using a Green function:
wherein u is s Is a scattering wave field, V is the scattering potential, u is the wave field value, G is the Green function,
is a position of the coordinate, and is,
to move the coordinate position, ω is the angular frequency.
Since the above formula is an expression of volume fraction, the scattering wave field is controlled by the scattering potential and the size of the scatterer.
In step 1, under the assumption of a first order bern Born approximation, i.e. weak scattering potential, the scatterers are small, so that the resulting scattered wavefield is much smaller than the background wavefield:
a scattered wave field under a first-order Born approximation is obtained:
wherein,
is a first-order scattered wave field, V is the scattering potential, u is the wave field value, G is the Green function,
to move the coordinate position, ω is the angular frequency,
and omega is an integral body for the coordinate position of the wave detection point.
The corresponding full wavefield is:
wherein u is the wave field value, u 0 Is the background wave field value, V is the scattering potential, G is the Green's function, ω is the angular frequency,
is the coordinate position of the wave detection point, omega is an integral body,
is a mobile coordinate location.
Substituting the formula (10) into the formula (8) to obtain a second-order scattered wave expression; repeating the above process to obtain a multi-order scattering wave field expression, and simply expressing the N-order scattering wave field by using a summation form:
wherein u is s In order to scatter the wave field,
for the ith order scattered wavefield, N is the total order,
is the coordinate position and ω is the angular frequency.
In step 2, when the scattering wave equation under the first-order Born approximation is adopted as a forward modeling operator, the LSRTM target functional under the L-2 norm is expressed as:
wherein J is a target functional,
in order to simulate the scattered wave field in the first order,
is a first-order scattered wave field, d n Is the nth order background wavefield.
The gradient of its corresponding velocity perturbation update is represented as:
wherein,
for speed perturbation, J is the target functional,
is the coordinate position, omega is the angular frequency,
in order to simulate the scattered wave field in the first order,
for first order simulation of the wave field, d n For 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 demodulator probe.
As can be seen from the analysis, the gradient includes two parts: simulating a gradient formed by a residual error between the background wave field and the simulated first-order scattering wave field and the first-order scattering wave field in the observation data:
wherein,
in order to have a first-order gradient,
is the coordinate position, omega is the angular frequency, u 0 Is the background wavefield value, G is the Green's function,
in order to detect the coordinate position of the point,
in order to move the coordinate position of the mobile phone,
in order to simulate the scattered wave field in the first order,
the wavefield is first order simulated.
And simulating a gradient formed by the background wave field and the observation data under the approximation of the background wave field and the high-order Born in the observation data:
wherein, g n Is the nth order gradient, ω is the angular frequency, u 0 In order to be the value of the background wave field,
is a position of the coordinate(s),
for the coordinate position of the detection point, d n For the nth order background wavefield, G is the Green function,
is a mobile coordinate location.
The first term of the gradient describes the velocity perturbation corresponding to the first-order scattered wave described by the forward operator, and 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 generate interference on the velocity perturbation gradient with the background wave field due to the fact that the other components cannot be accurately described.
In step 2, in order to suppress noise generated by the background wave field and the high-order scattered wave field in the observation data, a linear Berne approximation equation is added to the objective function of the nonlinear LSRTM as a constraint term, and the constraint of the first-order scattered wave in the nonlinear LSRTM is strengthened, so that a high-precision and high-resolution profile is obtained.
In step 2, the modified LSRTM objective function is represented as:
wherein, lambda is a constraint factor, the adjustment is carried out according to the intensity of a background wave field and a high-order scattered wave field in the observation data, J is a target functional,
is the coordinate position, omega is the angular frequency,
in order to simulate the scattered wave field in the first order,
for first order simulation of the wave field, d n For 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,
in order to detect the coordinate position of the point,
to move the coordinate position.
LSRTM after addition constraints is decomposed into two inverse problems: firstly, performing inversion on a first-order scattering wave field to obtain a relatively accurate first-order scattering wave field; then calculating the updating gradient of the speed disturbance on the basis; let the derivative of equation (16) for the first order scatter wavefield be 0, a reconstructed expression for the first order scatter wavefield can be obtained:
wherein, lambda is a constraint factor, the adjustment is carried out according to the intensity of a background wave field and a high-order scattered wave field in the observation data, J is a target functional,
is the coordinate position, omega is the angular frequency,
in order to simulate the scattered wave field in the first order,
for first order simulation of the wave field, d n For the nth order background wave field, u 0 Is the background wave field value, V is the scattering potential, G is the Green function,
in order to detect the coordinate position of the point,
to move the coordinate position.
Wherein, the lambda is a constraint factor,
is the coordinate position, omega is the angular frequency, u s 1 In order to simulate the scattered wave field in the first order,
for first order simulation of the wave field, d n For the nth order background wave field, u 0 Is the background wave field value, V is the scattering potential, G is the Green function,
is a mobile coordinate location.
As can be seen from the analysis (18), the reconstructed first-order scattered wave field is mainly composed of two parts: the proportion of the wave field formed by observation data and the simulated first-order scattering wave field is controlled by a constraint factor: when the constraint factor is larger, the influence of the first-order scattered wave on the reconstructed scattered wave field becomes larger, so that the influence of the background wave field and the high-order scattered wave on the first-order reconstructed scattered wave is weakened.
In step 2, after the scattering wave field is accurately reconstructed, the gradient of the target function to the velocity disturbance is obtained:
wherein, lambda is a constraint factor, J is a target functional,
is the coordinate position, omega is the angular frequency, u 0 Is the background wave field value, V is the scattering potential,
to move the coordinate position, G is the Green's function,
in order to be a speed disturbance,
for the first order of the simulated scattered wavefield, u is the wavefield value, u 0 Is the background wave field value.
And further updating the speed disturbance:
wherein alpha is the update step length, J is the target functional,
is a velocity disturbance. .
In step 3, for a general inverse problem, the convexity of the inverse problem is determined by the properties of the second derivative of the objective function pair and the inversion parameters, i.e. the Hessian matrix: when the Hessian matrix is determined, the Hessian matrix is semi-positive, namely the target function is fully convex, and the global extreme value of the target function is found by a local optimization method; conversely, the objective function is locally convex, and the inverse problem is affected by the local extremum and is difficult to converge to the global extremum.
In step 3, the Hessian matrix of the nonlinear LSRTM based on the linear berne approximation constraint on the first order scattered wavefield is:
wherein, lambda is a constraint factor, J is a target functional,
is a position of the coordinate, and is,
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 function,
in order to move the coordinate position of the mobile phone,
in order to simulate the scattered wave field in the first order,
for first order simulation of the wave field, d n For the n-th order background wave field,
is the coordinate position of the wave detection point.
The expression of the first-order scattering wave field shows that the wave field is influenced by model parameters, scattering potential and an observation system; wavefields under different observation systems are uncorrelated, then:
wherein H is a second order gradient, r is a position of a detection point, lambda is a constraint factor, J is a target functional,
is a position of the coordinate, and is,
is the coordinate position, omega is the angular frequency,
the scattered wavefield is modeled for the first order.
The nonlinear LSRTM based on the linear Bern approximate constraint shows strong convexity aiming at the first-order scattering wave field, and a more accurate first-order scattering wave field is obtained by a local optimization method; in addition, the first-order scattering wave field under the Born approximation is linear with the velocity disturbance, so that the accurate velocity disturbance updating gradient is obtained based on the reconstructed first-order scattering wave field.
The nonlinear least square reverse time migration method based on linear Bern approximation constraint in the invention proposes that linear Bern approximation is added into a target function of nonlinear LSRTM as a constraint term, and false images caused by weak scattering potential and small scatterer hypothesis are suppressed by enhancing the constraint of the linear Bern approximation in the inversion process. The linear berne approximation constraint-based inversion method decomposes the nonlinear LSRTM problem into two inversion problems: firstly, an optimal solution is found in a first-order scattering wave field space, namely, an accurate first-order scattering wave field is reconstructed as far as possible through observation data and a first-order scattering wave equation under the condition of given speed disturbance. An updated gradient of velocity perturbations is then calculated based on the first order scatter wavefield. The new objective function shows strong convexity on a Hessian matrix of the first-order scattering wave field, and the accuracy of the reconstruction of the first-order scattering wave field is ensured, so that the velocity disturbance updating gradient calculated according to the first-order scattering wave field expression is more accurate, and a high-precision and high-amplitude seismic profile can be obtained.
The invention has the beneficial effects that: the invention provides a nonlinear least square reverse time migration method based on linear Bern approximate constraint, which reconstructs an accurate first-order scattered wave field as much as possible through observation data and a first-order scattered wave equation. An updated gradient of velocity perturbations is then calculated based on the first-order scatter wavefield. The new objective function shows strong convexity on a Hessian matrix of the first-order scattering wave field, and the accuracy of the reconstruction of the first-order scattering wave field is ensured, so that the velocity disturbance updating gradient calculated according to the first-order scattering wave field expression is more accurate, and a high-precision and high-amplitude seismic profile can be obtained.
Drawings
FIG. 1 is a diagram of a Sigsbee2A model in an embodiment of the present invention;
FIG. 2 is a schematic illustration of inversion shot logging in accordance with an embodiment of the present invention 1;
fig. 3 is a schematic diagram of a sigbee 2A offset velocity model in embodiment 1 of the present invention;
FIG. 4 is a diagram of LSRTM results for 10 iterations in exemplary embodiment 1 of the present invention;
FIG. 5 is an enlarged view of the upper part of the salt of FIG. 4 in the embodiment 1 of the present invention;
FIG. 6 is an enlarged view of the lower portion of the salt of FIG. 4 according to embodiment 1 of the present invention;
FIG. 7 is a schematic illustration of inversion shot logging in accordance with an embodiment 2 of the present invention;
fig. 8 is a schematic diagram of a sigsebe 2A offset velocity model in embodiment 2 of the present invention;
FIG. 9 is a diagram of LSRTM results for 10 iterations in example 2 of the present invention;
FIG. 10 is an enlarged view of the upper part of the salt of FIG. 4 in the embodiment 2 of the present invention;
FIG. 11 is an enlarged view of the lower portion of the salt of FIG. 4 according to embodiment 2 of the present invention;
FIG. 12 is a schematic illustration of inversion shot records in example 3 of the present invention;
fig. 13 is a schematic diagram of a sigsebe 2A offset velocity model according to embodiment 3 of the present invention;
FIG. 14 is a diagram of LSRTM results for 10 iterations in example 3 of the present invention;
FIG. 15 is an enlarged view of the upper part of the salt of FIG. 4 in the embodiment 3 of the present invention;
FIG. 16 is an enlarged view of the lower portion of the salt of FIG. 4 according to embodiment 3 of the present invention;
FIG. 17 is a flowchart of an embodiment of the nonlinear least-squares reverse-time migration method based on the linear Bern's approximation constraint.
Detailed Description
It is to be understood that the following detailed description is exemplary and is intended to provide further explanation of the invention as claimed. Unless defined otherwise, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention belongs.
It is noted that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting of exemplary embodiments according to the invention. As used herein, the singular forms "a", "an" and "the" are intended to include the plural forms as well, and it should be understood that when the terms "comprises" and/or "comprising" are used in this specification, they specify the presence of the stated features, steps, operations, and/or combinations thereof, unless the context clearly indicates otherwise.
As shown in fig. 17, fig. 17 is a flowchart of the nonlinear least square reverse time migration method based on the linear berne approximation constraint of the present invention. The nonlinear least square reverse time migration method based on the linear Bern approximate constraint comprises the following steps:
101, under the condition of given speed disturbance, reconstructing an accurate first-order scattered wave field as far as possible through observation data and a first-order scattered wave equation;
the homogeneous Helmholtz equation for a scalar field of propagation in the subsurface medium is:
wherein,
is the coordinate position, omega is the angular frequency, delta is the hamiltonian, u is the wavefield value, Γ is the positive operator, v is the velocity value.
Introduction of background velocity field
We can get:
wherein,
is the coordinate position, omega is the angular frequency, delta is the Hamiltonian, u is the wave field value, v is the velocity value 0 Is the background velocity.
The scattering potential is defined as:
wherein,
is the coordinate position, omega is the angular frequency, V is the velocity value, V is the scattering potential, V 0 In the case of the background velocity,
i.e. the velocity disturbances of interest to the LSRTM. Depending on the scattering potential, the above equation can be simplified as:
wherein,
is the coordinate position, ω is the angular frequency, Δ is the Hamiltonian, u is the wave field value, v 0 For background velocity field, V is the scattering potential.
Assuming that the incident wavefield propagating in the background medium satisfies the homogeneous Helmholtz equation:
wherein,
is the coordinate position, omega is the angular frequency, delta is the Hamiltonian, u is the wave field value, v 0 Is the background velocity field.
Equation (5) is reduced to a non-homogeneous Helmholtz equation describing the scattered wavefield that results after the incident wavefield propagating in the background medium encounters the scatterers:
wherein,
is the coordinate position, omega is the angular frequency, delta is the Hamiltonian, u is the wave field value, v 0 As background velocity field, u s Is the scattering wavefield and V is the scattering potential.
Defining a Green function:
wherein,
is the coordinate position, omega is the angular frequency, delta is the Hamiltonian, G is the Green function, delta is the impulse function,
is a mobile coordinate location.
An analytical solution of the scattered wavefield may be obtained using the Green function:
wherein u is s Is a scattering wave field, V is the scattering potential, u is the wave field value, G is the Green function,
is a position of the coordinate, and is,
to move the coordinate position, ω is the angular frequency.
Since the above formula is an expression of volume fraction, the scattering wave field is controlled by the scattering potential and the size of the scatterer. Under the assumption of a first-order Born approximation, i.e. the scattering potential is weak, the scatterers are small, so that the resulting scattered wavefield is much smaller than the background wavefield:
we can obtain the scattered wavefield at the first-order Born approximation:
wherein,
is a first-order scattered wave field, V is the scattering potential, u is the wave field value, G is the Green function,
to move the coordinate position, ω is the angular frequency,
omega is an integral body for the coordinate position of the demodulator probe.
The corresponding full wavefield is:
wherein u is the wave field value, u 0 Is the background wave field value, V is the scattering potential, G is the Green's function, ω is the angular frequency,
is the coordinate position of the wave detection point, omega is an integral body,
is a mobile coordinate location.
By substituting the formula (10) into the formula (8), a second order scattered wave expression can be obtained. Repeating the above process to obtain a multi-order scattering wave field expression, and simply expressing the N-order scattering wave field by using a summation form:
wherein u is s In order to scatter the wave field,
for the ith order scattered wavefield, N is the total order,
is the coordinate position and ω is the angular frequency.
deriving an updated gradient of the computed velocity perturbation based on the first-order scattering wavefield, comprising:
when the scattered wave equation under the first-order Born approximation is taken as a forward operator, the LSRTM target functional under the L-2 norm can be expressed as:
wherein J is a target functional,
in order to simulate the scattered wave field in the first order,
is a first-order scattered wave field, d n Is the nth order background wavefield.
The gradient of its corresponding velocity perturbation update can be expressed as:
wherein,
for speed perturbation, J is the target functional,
is the coordinate position, omega is the angular frequency,
in order to simulate the scattered wave field in the first order,
for first order simulation of the wave field, d n For the nth order background wave field, u 0 Is the background wave field value, V is the scattering potential, G is the Green function,
is the coordinate position of the wave detection point.
The analysis shows that the gradient comprises two parts: and simulating a gradient formed by a residual error between the background wave field and the simulated first-order scattering wave field and the first-order scattering wave field in the observation data:
wherein,
in order to be a first-order gradient,
is the coordinate position, omega is the angular frequency, u 0 Is the background wavefield value, G is the Green's function,
in order to detect the coordinate position of the point,
in order to move the coordinate position of the mobile phone,
in order to simulate the scattered wave field in the first order,
the wavefield is first order simulated.
And simulating a gradient formed by the background wave field and the observation data under the approximation of the background wave field and the high-order Born in the observation data:
wherein, g n Is the nth order gradient, ω is the angular frequency, u 0 As a value of the background wave field,
is a position of the coordinate(s),
for the coordinate position of the detection point, d n For the nth order background wave field, G is the Green function,
is a mobile coordinate location.
It is obvious that 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, i.e. under the first-order Born approximation, other components of the observed data generate interference with the background wave field on the gradient of the velocity perturbation due to the fact that the other components cannot be accurately described. In order to suppress noise generated by a background wave field and a high-order scattered wave field in observation data, a linear Berne approximation equation is added into an objective function of the nonlinear LSRTM as a constraint term, and the constraint of a first-order scattered wave in the nonlinear LSRTM is strengthened to obtain a high-precision and high-resolution profile. The improved LSRTM objective function can be expressed as:
wherein, lambda is a constraint factor, the adjustment is carried out according to the intensity of a background wave field and a high-order scattered wave field in the observation data, J is a target functional,
is the coordinate position, omega is the angular frequency,
in order to simulate the scattered wave field in the first order,
for first order simulation of the wave field, d n For the nth order background wave field, u 0 Is the background wave field value, V is the scattering potential, G is the Green function,
in order to detect the coordinate position of the point,
to move the coordinate position.
Adding constrained LSRTM can be decomposed into two inverse problems: firstly, inverting a first-order scattered wave field to obtain a relatively accurate first-order scattered wave field; an updated gradient of the velocity disturbance is then calculated on this basis. Let equation (16) have a derivative of 0 for the first order scattering wavefield, a reconstructed expression for the first order scattering wavefield can be obtained:
wherein, lambda is a constraint factor, the adjustment is carried out according to the intensity of a background wave field and a high-order scattered wave field in the observation data, J is a target functional,
is the coordinate position, omega is the angular frequency,
in order to simulate the scattered wave field in the first order,
for first order simulation of the wave field, d n For the nth order background wave field, u 0 Is the background wave field value, V is the scattering potential, G is the Green function,
in order to detect the coordinate position of the point,
is a mobile coordinate location.
Wherein, the lambda is a constraint factor,
is the coordinate position, omega is the angular frequency,
in order to simulate the scattered wave field in the first order,
for first order simulation of the wave field, d n For the nth order background wave field, u 0 Is the background wave field value, V is the scattering potential, G is the Green function,
is a mobile coordinate location.
As can be seen from the analysis (18), the reconstructed first-order scattered wave field is mainly composed of two parts: the proportion of the wave field formed by observation data and the simulated first-order scattering wave field is controlled by a constraint factor: when the constraint factor is larger, the influence of the first-order scattered wave on the reconstructed scattered wave field becomes larger, so that the influence of the background wave field and the high-order scattered wave on the first-order reconstructed scattered wave is weakened.
After the scattering wavefield is accurately reconstructed, the gradient of the objective function to the velocity perturbation can be obtained:
wherein, lambda is a constraint factor, J is a target functional,
is the coordinate position, omega is the angular frequency, u 0 Is the background wave field value, V is the scattering potential,
for moving coordinate positions, G is the Green function,
in order to be a speed disturbance,
for the first order of the simulated scattered wavefield, u is the wavefield value, u 0 Is the background wave field value.
The speed disturbance can then be updated:
wherein alpha is the updating step length, J is the target functional,
is a velocity perturbation. .
103, discussing the convexity of an inverse problem of the nonlinear LSRTM based on the linear Bern approximate constraint; for a general inverse problem, the convexity of the inverse problem is generally determined by the nature of the second derivative (Hessian matrix) of the objective function pair and the inversion parameters: when the Hessian matrix is determined, the Hessian matrix is semi-positive, namely the target function is fully convex, and the global extreme value of the target function can be found by a local optimization method; conversely, the objective function is locally convex, and the inverse problem is affected by the local extremum and is difficult to converge to the global extremum. The Hessian matrix of the nonlinear LSRTM based on the linear berne approximation constraint on the first order scattered wavefield is:
wherein, lambda is a constraint factor, J is a target functional,
is a position of the coordinate(s),
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 function,
in order to move the coordinate position of the mobile terminal,
in order to simulate the scattered wave field in the first order,
for first order simulation of the wave field, d n For the n-th order background wave field,
is the coordinate position of the demodulator probe.
As can be seen from the expression of the first-order scattered wave field, the wave field is affected by the model parameters, the scattering potential and the observation system. Wavefields under different observation systems are uncorrelated, then:
wherein, lambda is a constraint factor, J is a target functional,
is a position of the coordinate, and is,
is the coordinate position, omega is the angular frequency,
the scattered wavefield is modeled for the first order.
Namely, the nonlinear LSRTM based on the linear Berne approximate constraint shows strong convexity aiming at the first-order scattered wave field, and a more accurate first-order scattered wave field can be obtained by using a local optimization method. Furthermore, the first-order scattering wave field under the Born approximation is linear with the velocity perturbation, so we can get accurate velocity perturbation update gradients based on the reconstructed first-order scattering wave field.
And step 104, updating the gradient according to the velocity disturbance calculated by the first-order scattering wave field expression to obtain a seismic section with high precision and high amplitude preservation.
Example 1:
in the specific embodiment 1 to which the present invention is applied, the present invention is based on the least square reverse time migration of the primary scattered wave constraint, and is firstly used for testing the adaptability of the Sigsbee2A model (FIG. 1). The Sigsbee2A model contains a high-speed salt dome, so that the assumption of weak scattering potential and small scatterers cannot be met, and the effectiveness of the method can be well illustrated. The model size is 401 × 1601, the grid spacing is 21m longitudinally, and 22m transversely. And synthesizing the seismic record by adopting a time domain sound wave equation, wherein the sampling time is 7s, and the sampling interval is 1ms. 80 cannons are placed on the surface, the shot spacing is 440m, and the surface is received in a full aperture mode. Fig. 2 shows a shot record for inversion, the composite record contains multiple wave phenomena such as direct wave, rotating wave, etc., and the shot record also contains multiple waves such as interbed multiples, higher-order scattered waves, etc. due to the existence of salt dome and the complexity of the deep layer model. The background velocities used for the inversion are shown in fig. 3, with the model having no significant scatterers except for the large salt dome. During inversion testing, in order to reflect the importance of the constraint factor on inversion, common LSRTM, LSRTM with constraint factor of 1000 and LSRTM with constraint factor of 100000 are respectively used for inversion, the final inversion result is shown in FIG. 4, and FIG. 4 (a) is a conventional LSRTM result; (b) LSRTM results based on first order scattered wave constraints (constraint factor 1000); (c) Is LSRTM result based on first order scattered wave constraint (constraint factor 100000). Also, the inversion result of the deep layer is not obvious due to the over-strong energy of the direct wave and the rotating wave.
In order to compare the imaging results of the two methods more clearly. First, a portion of the salt was cut for comparison: due to the influence of direct waves and rotating waves, a large amount of interference noise exists in an inversion result given by the common LSRTM, and the shallow layer arc drawing phenomenon is serious. After the constraint of the primary scattered wave is introduced, the interference noise is obviously weakened, and the resolution of the inversion result is further improved. With the further enhancement of the constraint of the first-order scattered waves, noises generated by background wave fields such as direct waves, rotating waves and the like are basically eliminated, and the shallow layer position of the inversion result is more definite. Fig. 5 is an enlarged view of the upper portion of the salt for fig. 4.
Then, the next part of the method for salt is selected for accurate comparison: due to the complexity of observation data, a common LSRTM result contains a large amount of noise and does not invert scattering points, the inversion result of the LSRTM with the added constraint has better fidelity and improved resolution, but does not invert scattering points; by further enforcing the constraints of the first-order scattered waves, LSRTM ultimately gives better inversion results: the position of the horizon and scattering points under the salt is accurate, and the resolution is high. Fig. 6 is an enlarged view of the lower portion of the salt for fig. 4.
Least squares reverse time migration can give high resolution and high amplitude preservation of imaging results by minimizing the difference between the simulated data and the observed data, and is less prone to low frequency noise and artifacts than conventional reverse time migration methods. The LSRTM theory is based on the assumption of first-order scattered waves under Born approximation, and its corresponding scattered wave expression can only describe weak scattering potentials and wave fields under small scatterers. However, the actual propagation process of the seismic waves is a nonlinear process and comprises various wave phenomena, the underground medium is very complex, and the conditions of weak scattering potential and small scattering bodies are generally difficult to meet. When full wave field information is used for inversion, high-order scattered waves generated by strong scattering potential and large scatterers, direct waves, rotating waves and other waves generate false images with first-order scattered waves, and imaging quality is affected.
Therefore, in view of the drawbacks of the above-mentioned techniques, the present invention proposes to add the expression of the linear berne approximation as a constraint term to the objective function of the nonlinear LSRTM, suppressing artifacts caused by weak scattering potentials and small scatterer hypotheses by enhancing the constraint of the first-order scattered waves during the inversion process.
Example 2:
in the embodiment 2 to which the present invention is applied, the present invention is based on the least square reverse time migration of the primary scattered wave constraint, and is firstly used for testing the adaptability of the Sigsbee2A model (figure 1). The Sigsbee2A model contains a high-speed salt dome, so that the assumption of weak scattering potential and small scatterers cannot be met, and the effectiveness of the method can be well illustrated. The model size was 201 × 801, the grid spacing was 42m in the vertical direction and 44m in the horizontal direction. And synthesizing the seismic record by adopting a time domain sound wave equation, wherein the sampling time is 7s, and the sampling interval is 1ms. 80 cannons are placed on the surface, the shot spacing is 440m, and the surface is received in a full aperture mode. Fig. 7 shows a shot record for inversion, which contains multiple wave phenomena such as direct waves, gyrotrons, etc., and also includes multiple wave phenomena such as interbed multiples, higher-order scattered waves, etc., due to the presence of salt domes and the complexity of deep models. The background velocities used for the inversion are shown in fig. 8, and the model has no significant scatterers except for large salt domes. During inversion testing, in order to reflect the importance of the constraint factor to inversion, common LSRTM, LSRTM with constraint factor of 1000 and LSRTM with constraint factor of 100000 are respectively used for inversion, and the final inversion result is shown in fig. 9, where fig. 9 (a) is a conventional LSRTM result; (b) LSRTM results based on first order scattered wave constraints (constraint factor 1000); (c) Is LSRTM result based on first order scattered wave constraint (constraint factor 100000). Also, the inversion result of the deep layer is not obvious due to the over-strong energy of the direct wave and the rotating wave.
In order to compare the imaging results of the two methods more clearly. First, a portion of the salt was cut off for comparison: due to the influence of direct waves and rotating waves, a large amount of interference noise exists in an inversion result given by the common LSRTM, and the shallow arc drawing phenomenon is serious. After the constraint of the primary scattered wave is introduced, the interference noise is obviously weakened, and the resolution of the inversion result is further improved. With the constraint of the first-order scattered waves strengthened further, noise generated by background wave fields such as direct waves, rotating waves and the like is basically eliminated, and the shallow layer position of an inversion result is more definite. Fig. 12 is an enlarged view of the upper portion of the salt for fig. 9.
Then, the next part of the method for salt is selected for accurate comparison: due to the complexity of observation data, a common LSRTM result contains a large amount of noise and does not invert scattering points, the inversion result of the LSRTM with the added constraint has better fidelity and improved resolution, but does not invert scattering points; by further enforcing the constraints of the first-order scattered waves, LSRTM ultimately gives better inversion results: the positions of the horizon and scattering points under the salt are accurate, and the resolution ratio is high. Fig. 11 is an enlarged view of the lower portion of the salt for fig. 9.
Example 3:
in the specific embodiment 3 to which the present invention is applied, the present invention is based on the least square reverse time migration of the primary scattered wave constraint, and is firstly used for testing the adaptability of the Sigsbee2A model (figure 1). The Sigsbee2A model contains a high-speed salt dome, so that the assumption of weak scattering potential and small scatterers cannot be met, and the effectiveness of the method can be well illustrated. The model size was 101 × 401, the grid spacing was 84m in the vertical direction and 88m in the horizontal direction. And synthesizing the seismic record by adopting a time domain sound wave equation, wherein the sampling time is 7s, and the sampling interval is 1ms. 80 cannons are placed on the surface, the interval of the cannon points is 440m, and the surface is received in a full aperture mode. Fig. 8 shows a shot record for inversion, the composite record contains multiple wave phenomena such as direct wave, rotating wave, etc., and the shot record also contains multiple waves such as interbed multiples, higher-order scattered waves, etc. due to the existence of salt dome and the complexity of the deep model. The background velocities used for the inversion are shown in fig. 14, with the model having no significant scatterers except for the large salt dome. During inversion testing, in order to reflect the importance of the constraint factor on inversion, common LSRTM, LSRTM with constraint factor of 1000 and LSRTM with constraint factor of 100000 are respectively used for inversion, the final inversion result is shown in fig. 14, and fig. 14 (a) is a conventional LSRTM result; (b) LSRTM results based on first order scattered wave constraints (constraint factor 1000); (c) Is LSRTM result based on first order scattered wave constraint (constraint factor 100000). Also, the inversion result of the deep layer is not obvious due to the over-strong energy of the direct wave and the rotating wave.
In order to compare the imaging results of the two methods more clearly. First, a portion of the salt was cut off for comparison: due to the influence of direct waves and rotating waves, a large amount of interference noise exists in an inversion result given by the common LSRTM, and the shallow layer arc drawing phenomenon is serious. After the constraint of the primary scattered wave is introduced, the interference noise is obviously weakened, and the resolution of the inversion result is further improved. With the constraint of the first-order scattered waves strengthened further, noise generated by background wave fields such as direct waves, rotating waves and the like is basically eliminated, and the shallow layer position of an inversion result is more definite. Fig. 15 is an enlarged view of the upper portion of the salt from fig. 14.
Then, the next part of the salt method is selected for accurate comparison: due to the complexity of observation data, a common LSRTM result contains a large amount of noise and does not invert scattering points, while the inversion result of the LSRTM with increased constraints has better fidelity and improved resolution, but does not invert scattering points; by further enforcing the constraints of the first-order scattered waves, LSRTM ultimately gives a good inversion result: the positions of the horizon and scattering points under the salt are accurate, and the resolution ratio is high. Fig. 16 is an enlarged view of the lower portion of the salt from fig. 14.
Finally, it should be noted that: although the present invention has been described in detail with reference to the foregoing embodiments, it will be apparent to those skilled in the art that changes may be made in the embodiments and/or equivalents thereof without departing from the spirit and scope of the invention. Any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.
Except for the technical features described in the specification, the method is known by the technical personnel in the field.
Claims (10)
1. The nonlinear least square reverse time migration method based on the linear Berne approximate constraint is characterized by comprising the following steps of:
step 1, reconstructing a first-order scattered wave field through observation data and a first-order scattered wave equation under the condition of given speed disturbance;
step 2, deducing an updating gradient of the calculated speed disturbance based on the first-order scattering wave field;
step 3, discussing the convexity of an inverse problem of the nonlinear LSRTM based on the linear Berne approximate constraint;
and 4, updating the gradient according to the velocity disturbance calculated by the first-order scattering wave field expression to obtain a seismic section with high precision and high amplitude preservation.
2. The nonlinear least squares reverse time migration method based on linear berne approximation constraint as claimed in claim 1, wherein in step 1, the homogeneous Helmholtz equation corresponding to the propagated scalar field in the subsurface medium is:
wherein,
is a coordinate position, omega is an angular frequency, delta is a Hamiltonian, u is a wave field value, gamma is a positive operator, and v is a velocity value;
introduction of background velocity field
Obtaining:
wherein,
is the coordinate position, omega is the angular frequency, delta is the Hamiltonian, u is the wavefield value, v is the velocity value 0 As background velocity;
the scattering potential is defined as:
wherein,
is the coordinate position, omega is the angular frequency, V is the velocity value, V is the scattering potential, V 0 In the case of the background velocity,
namely the speed disturbance interested by LSRTM;
according to the scattering potential, the above formula is simplified as follows:
wherein,
is the coordinate position, ω is the angular frequency, Δ is the Hamiltonian, u is the wave field value, v 0 For background velocity field, V is the scattering potential.
3. The nonlinear least squares reverse time migration method based on linear berne approximation constraint according to claim 2, characterized in that in step 1, the incident wavefield propagating in the background medium is assumed to satisfy the homogeneous Helmholtz equation:
wherein,
is the coordinate position, omega is the angular frequency, delta is the Hamiltonian, u is the wave field value, v 0 Is the background velocity field;
equation (5) is reduced to a non-homogeneous Helmholtz equation describing the scattered wavefield that results after the incident wavefield propagating in the background medium encounters the scatterers:
wherein,
is the coordinate position, omega is the angular frequency, delta is the Hamiltonian, u is the wave field value, v 0 As background velocity field, u s Is a scattering wave field, V is a scattering potential;
defining a Green function:
wherein,
is the coordinate position, omega is the angular frequency, delta is the Hamiltonian, G is the Green function, delta is the impulse function,
is a mobile coordinate position;
obtaining an analytic solution of the scattering wave field by using a Green function:
wherein u is s Is a scattering wave field, V is the scattering potential, u is the wave field value, G is the Green function,
is a position of the coordinate, and is,
is the mobile coordinate position, omega is the angular frequency;
since the above formula is an expression of volume fraction, the scattering 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 Born approximation constraint according to claim 3, characterized in that in step 1, under the assumption of first-order Born approximation, i.e. scattering potential is weak, scatterers are small, so that the generated scattering wave field is much smaller than the background wave field:
a scattered wavefield at a first-order Born approximation is obtained:
wherein,
is a first-order scattered wave field, V is the scattering potential, u is the wave field value, G is the Green function,
to move the coordinate position, ω is the angular frequency,
is the coordinate position of the demodulator probe, and omega is an integral body;
the corresponding full wavefield is:
wherein u is the wave field value, u 0 Is the background wave field value, V is the scattering potential, G is the Green's function, ω is the angular frequency,
is the coordinate position of the wave detection point, omega is an integral body,
is a mobile coordinate position;
substituting the formula (10) into the formula (8) to obtain a second-order scattered wave expression; repeating the above process to obtain a multi-order scattering wave field expression, and simply expressing the N-order scattering wave field by using a summation form:
wherein u is s In order to scatter the wave field,
is the ith order scattered wave field, N is the total order,
is the coordinate position and ω is the angular frequency.
5. The nonlinear least squares reverse time migration method based on linear Born approximation constraint of claim 1, characterized in that in step 2, when the scattered wave equation under first order Born approximation is adopted as the forward modeling operator, the LSRTM target functional under L-2 norm is expressed as:
wherein J is a target functional,
in order to simulate the scattered wave field in the first order,
is a first-order scattered wave field, d n Is the nth order background wave field;
the gradient of its corresponding velocity perturbation update is represented as:
wherein,
for the velocity perturbation, J is the target functional,
is the coordinate position, omega is the angular frequency,
in order to simulate the scattered wave field in the first order,
for first order simulation of the wave field, d n For the nth order background wave field, u 0 Is the background wave field value, V is the scattering potential, G is the Green function,
is the coordinate position of the wave detection point;
the analysis shows that the gradient comprises two parts: and simulating a gradient formed by a residual error between the background wave field and the simulated first-order scattering wave field and the first-order scattering wave field in the observation data:
wherein,
in order to be a first-order gradient,
is the coordinate position, omega is the angular frequency, u 0 Is the background wave field value, G is the Green function,
in order to detect the coordinate position of the point,
in order to move the coordinate position of the mobile terminal,
in order to simulate the scattered wave field in the first order,
simulating a wavefield for a first order;
and simulating a gradient formed by the background wave field and the observation data under the approximation of the background wave field and the high-order Born in the observation data:
wherein, g n Is the nth order gradient, ω is the angular frequency, u 0 In order to be the value of the background wave field,
is a position of the coordinate, and is,
for the coordinate position of the detection point, d n For the nth order background wavefield, G is the Green function,
is a mobile coordinate position;
the first term of the gradient describes the velocity perturbation corresponding to the first-order scattered wave described by the forward operator, and 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 generate interference on the velocity perturbation gradient with the background wave field due to the fact that the other components cannot be accurately described.
6. The nonlinear least-squares reverse-time migration method based on linear born approximation constraint of claim 5, wherein 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 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 the high-precision and high-resolution profile.
7. The nonlinear least squares reverse time migration method based on linear berne approximation constraint according to claim 6, characterized in that in step 2, the modified LSRTM objective function is expressed as:
wherein, lambda is a constraint factor, the adjustment is carried out according to the intensity of a background wave field and a high-order scattered wave field in the observation data, J is a target functional,
is the coordinate position, omega is the angular frequency,
in order to simulate the scattered wave field in the first order,
for first order simulation of the wave field, d n For the nth order background wave field, u 0 Is the background wave field value, V is the scattering potential, G is the Green function,
in order to detect the coordinate position of the point,
is a mobile coordinate position;
LSRTM after addition constraints is decomposed into two inverse problems: firstly, performing inversion on a first-order scattering wave field to obtain a relatively accurate first-order scattering wave field; then calculating the updating gradient of the speed disturbance on the basis; let the derivative of equation (16) for the first order scatter wavefield be 0, a reconstructed expression for the first order scatter wavefield can be obtained:
wherein, lambda is a constraint factor, the adjustment is carried out according to the intensity of a background wave field and a high-order scattering wave field in the observation data, J is a target functional,
is the coordinate position, omega is the angular frequency,
in order to simulate the scattered wave field in the first order,
for first order simulation of the wave field, d n For the nth order background wave field, u 0 Is the background wave field value, V is the scattering potential, G is the Green function,
in order to detect the coordinate position of the point,
is a mobile coordinate position;
wherein, the lambda is a constraint factor,
is the coordinate position, omega is the angular frequency,
in order to simulate the scattered wave field in the first order,
for first order simulation of the wave field, d n For 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 a mobile coordinate position;
as can be seen from the analysis (18), the reconstructed first-order scattered wave field is mainly composed of two parts: the proportion of the wave field formed by observation data and the simulated first-order scattering wave field is controlled by a constraint factor: when the constraint factor is larger, the influence of the first-order scattered waves on the reconstructed scattered wave field becomes larger, so that the influence of the background wave field and the high-order scattered waves on the first-order reconstructed scattered waves is weakened.
8. The nonlinear least squares reverse time migration method based on linear berne approximation constraint according to claim 7, characterized in that in step 2, after the accurate reconstruction of the scattered wave field, the gradient of the objective function to velocity perturbation is obtained:
wherein, lambda is a constraint factor, J is a target functional,
is the coordinate position, omega is the angular frequency, u 0 Is the background wave field value, V is the scattering potential,
to move the coordinate position, G is the Green's function,
in order to be a speed disturbance,
for first order simulation of the scattered wavefield, u is the wavefield value 0 Is the background wave field value;
and further updating the speed disturbance:
wherein alpha is the updating step length, J is the target functional,
is a velocity perturbation.
9. The nonlinear least squares reverse time migration method based on linear berne approximation constraint according to claim 1, wherein in step 3, for a general inverse problem, the property of the Hessian matrix, which is the second derivative of the objective function pair and the inversion parameter, determines the convexity of the inverse problem: when the Hessian matrix is determined, the Hessian matrix is semi-positive, namely the target function is fully convex, and the global extremum of the target function is found by a local optimization method; on the contrary, the objective function is locally convex, and the inverse problem is affected by the local extremum and is difficult to converge to the global extremum.
10. The nonlinear least squares reverse time migration method based on linear born approximation constraint of claim 9, wherein in step 3, the Hessian matrix of nonlinear LSRTM based on linear born approximation constraint on first order scattered wavefields is:
wherein, lambda is a constraint factor, J is a target functional,
is a position of the coordinate, and is,
is the coordinate position, omega is the angular frequency, u 0 Is the background wave field value, V is the scattering potential, G is the Green function,
in order to move the coordinate position of the mobile terminal,
in order to simulate the scattered wave field in the first order,
for first order simulation of the wave field, d n For the n-th order background wave field,
is the coordinate position of the wave detection point;
the expression of the first-order scattering wave field shows that the wave field is influenced by model parameters, scattering potential and an observation system; wavefields under different observation systems are uncorrelated, then:
wherein, lambda is a constraint factor, J is a target functional,
is a position of the coordinate, and is,
is the coordinate position, omega is the angular frequency,
simulating a scattered wave field for a first order;
the nonlinear LSRTM based on the linear Bern approximate constraint shows strong convexity aiming at the first-order scattering wave field, and a more accurate first-order scattering wave field is obtained by a local optimization method; furthermore, the first-order scattering wave field under the Born approximation is linear with the velocity perturbation, so that an accurate velocity perturbation update gradient is obtained based on the reconstructed first-order scattering wave field.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110602166.0A CN115480294A (en) | 2021-05-31 | 2021-05-31 | Nonlinear least square reverse time migration method based on linear Beren approximate constraint |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110602166.0A CN115480294A (en) | 2021-05-31 | 2021-05-31 | Nonlinear least square reverse time migration method based on linear Beren approximate constraint |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN115480294A true CN115480294A (en) | 2022-12-16 |
Family
ID=84419514
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110602166.0A Pending CN115480294A (en) | 2021-05-31 | 2021-05-31 | Nonlinear least square reverse time migration method based on linear Beren approximate constraint |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115480294A (en) |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20170115418A1 (en) * | 2014-06-30 | 2017-04-27 | Cgg Services Sas | Seismic data processing using matching filter based cost function optimization |
| CN107783190A (en) * | 2017-10-18 | 2018-03-09 | 中国石油大学(北京) | A kind of least square reverse-time migration gradient updating method |
| KR20210042720A (en) * | 2019-10-10 | 2021-04-20 | 인하대학교 산학협력단 | Seismic waveform inversion method and apparatus for underground velocity model |
-
2021
- 2021-05-31 CN CN202110602166.0A patent/CN115480294A/en active Pending
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20170115418A1 (en) * | 2014-06-30 | 2017-04-27 | Cgg Services Sas | Seismic data processing using matching filter based cost function optimization |
| CN107783190A (en) * | 2017-10-18 | 2018-03-09 | 中国石油大学(北京) | A kind of least square reverse-time migration gradient updating method |
| KR20210042720A (en) * | 2019-10-10 | 2021-04-20 | 인하대학교 산학협력단 | Seismic waveform inversion method and apparatus for underground velocity model |
Non-Patent Citations (9)
| Title |
|---|
| KE CHEN 等: "Time-domain elastic Gauss–Newton full-waveform inversion: a matrix-free approach", GEOPHYSICAL JOURNAL INTERNATIONAL, vol. 223, no. 2, 21 July 2020 (2020-07-21), pages 1007 - 1039 * |
| 任志明: "声波和弹性波波动方程有限差分正反演方法研究", 中国博士学位论文全文数据库 基础科学辑, no. 2, 15 February 2018 (2018-02-15), pages 011 - 80 * |
| 刘玉金 等: "扩展成像条件下的最小二乘逆时偏移", 地球物理学报, vol. 58, no. 10, 31 October 2015 (2015-10-31), pages 3771 - 3782 * |
| 刘玉金: "复杂介质条件下基于反演的地震成像方法研究", 中国博士学位论文全文数据库 基础科学辑, no. 1, 15 January 2018 (2018-01-15), pages 011 - 129 * |
| 孙小东 等: "基于共轭梯度法和互相关的最小二乘逆时偏移及应用(英文)", 应用地球物理(英文版), vol. 14, no. 003, 30 September 2017 (2017-09-30), pages 381 - 386 * |
| 徐凯: "基于 GPU/CPU 协同加速的最小二乘逆时偏移", 中国优秀硕士学位论文全文数据库 基础科学辑, no. 4, 15 April 2018 (2018-04-15), pages 011 - 563 * |
| 方修政 等: "基于逆散射成像条件的最小二乘逆时偏移", 地球物理学报, vol. 61, no. 09, 30 September 2018 (2018-09-30), pages 3770 - 3782 * |
| 李武群 等: "二阶Born近似及其在逆散射GRT反演中的应用", 地球物理学进展, vol. 32, no. 02, 30 April 2017 (2017-04-30), pages 645 - 656 * |
| 王华忠 等: "最小二乘叠前深度偏移成像理论与方法", 石油物探, vol. 56, no. 02, 31 March 2017 (2017-03-31), pages 159 - 170 * |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US10324205B2 (en) | Device and method for full waveform inversion | |
| AU2010292176B2 (en) | Dip guided full waveform inversion | |
| AU2007302695B2 (en) | Iterative inversion of data from simultaneous geophysical sources | |
| CA2972033C (en) | Multistage full wavefield inversion process that generates a multiple free data set | |
| CN110187382B (en) | Traveling time inversion method for wave equation of reverse wave and reflected wave | |
| US20130258810A1 (en) | Method and System for Tomographic Inversion | |
| RU2570827C2 (en) | Hybrid method for full-waveform inversion using simultaneous and sequential source method | |
| CN110780341B (en) | Anisotropic seismic imaging method | |
| GB2584196A (en) | Methods and devices performing adaptive quadratic Wasserstein full-waveform inversion | |
| CN107076868A (en) | Declined using many subwaves of imperfect geological data and drop and be imaged | |
| CN115480294A (en) | Nonlinear least square reverse time migration method based on linear Beren approximate constraint | |
| Hadden et al. | Anisotropic full-waveform inversion of crosshole seismic data: A vertical symmetry axis field data application | |
| CN111025393B (en) | Reservoir prediction method, device, equipment and medium for stratum containing thin coal seam | |
| CN111175822B (en) | Strong scattering medium inversion method for improving direct envelope inversion and disturbance decomposition | |
| Lecomte | Hybrid modeling with ray tracing and finite difference | |
| CN111665549A (en) | Inversion method of stratum acoustic wave attenuation factor | |
| CN111665550A (en) | Underground medium density information inversion method | |
| CN112462428B (en) | Multi-component seismic data migration imaging method and system | |
| Galuzzi | Modelling and optimization techniques for acoustic Full Waveform Inversion in seismic exploration | |
| CN118534532A (en) | Method, equipment and storage medium for modeling inversion of elastic emission waveforms of seismic exploration | |
| Chen | Improving the Quality of Seismic Images by Deterministic Inversion and Machine Learning Methods | |
| CN115685322A (en) | Viscous medium target area imaging method based on virtual detection point dyeing algorithm | |
| CN111665552A (en) | Acoustic parameter acquisition method for landslide risk evaluation | |
| CN111665551A (en) | Acoustic parameter acquisition method for bridge foundation detection | |
| Li | Data Assimilation Methods in Seismic Inversion |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |