KR101820850B1

KR101820850B1 - Seismic imaging apparatus and method using iterative direct waveform inversion

Info

Publication number: KR101820850B1
Application number: KR1020150102111A
Authority: KR
Inventors: 신창수
Original assignee: 서울대학교 산학협력단
Priority date: 2014-07-25
Filing date: 2015-07-20
Publication date: 2018-01-22
Also published as: KR20160012922A

Abstract

The present invention relates to an acoustic wave imaging technique and relates to an imaging technique for modeling an underground structure by updating a velocity model of each frequency band in ascending order of frequency, descending order of frequency, or arbitrary order.
One challenge that the present invention addresses is to directly calculate the difference between the actual surface velocity and the initial estimated velocity.
In order to visualize the subterranean structure of a measurement area according to one aspect, an elastic wave imaging device using repetitive application of direct waveform inversion utilizes a virtual scattering source and a parameter perturbation calculated from an updated reference wave field A waveform inversion section for updating the reference velocity model by updating the reference velocity model updated in the previous frequency band in the ascending order of frequency, descending order of frequency, or arbitrary order to the reference velocity model of the next frequency band, And an underground structure imaging unit for imaging the underground structure from the model.

Description

[0001] SEISMIC IMAGING APPARATUS AND METHOD FOR USING ITERATED DIRECT WAVEFORM INVERSION [0002]

The present invention relates to an elastic wave imaging technique, and relates to an imaging technique for modeling an underground structure by updating a velocity model of each frequency band while shifting the frequency band.

Seismic tomography is a means of imaging the interior of the earth. Geophysicists have been studying seismic tomography to investigate the global scale structure within the earth through a passive method of analyzing signals from earthquakes. Exploration geophysicists have used seismic tomography to look for gas or oil reservoirs from the top of the earth from seismic exploration data.

Tomography methods are commonly used to create an initial velocity model for migration in geophysical exploration. In many available seismic tomography methods, waveform tomography uses seismic data as much as possible. Waveform tomography is known as full waveform inversion, and has continued to evolve in computing technology.

Full waveform inversion is a method of repeatedly calculating the ground velocity structure using seismic data. In order to find the solution by the complete waveform inversion, it is necessary to change the velocity model repeatedly to minimize the objective function defined by the difference between the acquired data and the composite data, and many methodologies have been proposed to implement it (Menke, 1984; Tarantola, 2005).

However, the complete waveform inversion based on generalized inversion theory requires a large computational computation time and storage space for jacobian matrix computation, hessian matrix composition, inverse matrix or determinant solution for multivariable functions. It has limitations.

One challenge that the present invention addresses is to directly calculate the difference between the actual surface velocity and the initial estimated velocity.

Another problem to be solved by the present invention is that the steepest descent method is not used.

Another problem to be solved by the present invention is to estimate the difference in the physical property between the initial velocity model and the inverse object medium by using the Lippmann-Schwinger equation.

In order to visualize the subterranean structure of a measurement area according to one aspect, an elastic wave imaging device using repetitive application of direct waveform inversion utilizes a virtual scattering source and a parameter perturbation calculated from an updated reference wave field A waveform inversion section for updating the reference velocity model by shifting the frequency band in the set order and updating the reference velocity model updated in the previous frequency band as the reference velocity model for the next frequency band, Lt; RTI ID = 0.0 > imaging < / RTI >

In another aspect, the waveform inversion section sets the updated reference speed model in the last frequency band back to the reference speed model in the first frequency band, sets the updated speed model in the previous frequency band as the reference speed model in the next frequency band .

In another aspect, if the frequency is a complex value, the waveform inversion section can shift the frequency band while changing the damping constant for adjusting the imaginary part and the magnitude of the real part in the set order.

In another aspect, the set order is an ascending order, a descending order, or an arbitrary order.

In another aspect, the waveform inversion section may include a reference wave length calculating section for calculating a reference wave length, which is a solution of the wave equation, using data collected from at least one of a transmitter and a receiver, A scattered wave field calculating unit for calculating a scattered wave field from the virtual scattering transmission source, a reference wave field updating unit for updating the reference wave field by adding the reference wave field and the scattered wave field, A parameter perturbation calculation unit for calculating a parameter perturbation from the virtual scattering transmission source and the updated reference wave field, a reference speed model updating unit for updating the reference speed model from the data collected from at least one of the transmitter and the receiver, and the calculated parameter perturbation .

In another aspect, the reference wave length calculation unit may calculate a reference wave length field by using a finite element method or a finite difference method in a frequency domain waveform equation transformed from a time domain waveform equation through a Fourier transform have.

In another aspect, the virtual scatter transmission source calculation unit may calculate a virtual scatter transmission source using a function of residual data and Green in an expression derived from a Lippmann Schwinger equation.

In yet another aspect, the parameter perturbation calculator may calculate the parameter perturbation by applying a Newton method or a least squares method to the objective function.

According to the proposed invention, since the steepest descent method which is repeatedly calculated to reduce the residual data is not used, the calculation cost required to perform the waveform inversion can be reduced.

Furthermore, according to the proposed invention, the step length for updating the velocity model by directly inversely interpolating the parameter perturbation may not be used.

Furthermore, by using the sequential frequency strategy, a precise speed model with small or large differences in the material properties between the initial velocity model and the inverse target medium can be obtained.

Further, according to the proposed invention, the updated speed model is shifted to the initial speed model by shifting the frequency band, and the inverse calculation is performed to obtain the final speed model close to the actual medium.

FIG. 1 shows a schematic configuration of an elastic-wave imaging apparatus using repetitive application of direct waveform inversion according to an embodiment.
Figure 2 shows the principle of waveform inversion according to one embodiment.
FIG. 3 shows a detailed configuration of an elastic-wave imaging apparatus using repetitive application of direct waveform inversion according to an embodiment.
4 is a flowchart of an elastic wave imaging method using repetitive application of direct waveform inversion according to an embodiment.
5 is a detailed flowchart of the waveform inversion step according to one embodiment.
Figure 6 shows the Marmousi velocity model.
7 (a) shows a shot-gather seismogram and (b) shows a frequency spectrum of a Marmousi data set.
8 shows an initial velocity model.
Figures 9 (a), 9 (b) and 9 (c) each show the reference wave field propagated at 4.8 Km shot position in the 5.0 Hz, 10.0 Hz and 15.0 Hz initial velocity models.
Each of Figures 10 (a), (b), and (c) shows a virtual scattering source obtained from residual data corresponding to a point source at a shot position of 4.8 Km at 5.0 Hz, 10.0 Hz and 15.0 Hz.
Each of Figures 11 (a), (b) and (c) shows the scattered wave field generated by the virtual scattering source shown in Figure 10 at 5.0 Hz, 10.0 Hz and 15.0 Hz.
Each of Figs. 12 (a), (b) and (c) shows the generated wave at the shot position of 4.8 Km, which is the sum of the reference wave field of Fig. 9 and the scattered wave field of Fig. 10 at 5.0 Hz, 10.0 Hz and 15.0 Hz Field.
Figures 13 (a), 13 (b) and 13 (c) each show the actual wave field at a shot position of 4.8 Km at 5.0 Hz, 10.0 Hz and 15.0 Hz.
Figures 14 (a), 14 (b) and 14 (c) show parameter perturbations at 5.0 Hz, 10.0 Hz and 15.0 Hz, respectively.
Figure 15 shows the inversed velocity model.
16 (a), 16 (b), and 16 (c) are diagrams each showing the Shotgarder seismic data of synthesized field data, the synthesized seismic data in the initial velocity model, and the seismic data obtained by direct waveform inversion / RTI >
17 (a) and 17 (b) show Shotgarder seismic data of a noise-mixed Marmousi data set with signal-to-noise ratios of 12.75 and 25.50, respectively.
18 (a) and 18 (b) show the inverse model of the noise in the noise-mixed Marmousi data set with S / N ratios of 12.75 and 25.50, respectively.

The foregoing and further aspects are embodied through the embodiments described with reference to the accompanying drawings. It is to be understood that the components of each embodiment are capable of various combinations within an embodiment as long as no other mention or mutual contradiction exists. Furthermore, the proposed invention may be embodied in many different forms and is not limited to the embodiments described herein.

In order to clearly illustrate the claimed invention, parts not related to the description are omitted, and like reference numerals are used for like parts throughout the specification. And, when a section is referred to as "including " an element, it does not exclude other elements unless specifically stated to the contrary. As used herein, the term " block " refers to a block of hardware or software configured to be changed or pluggable, i.e., a unit or block that performs a specific function in hardware or software.

FIG. 1 shows a schematic configuration of an underground structure imaging apparatus according to an embodiment of the present invention.

1, an underground structure imaging apparatus may include a transmission source 102, a receiver 103, and a signal processing apparatus 104 installed in a measurement target area 101. FIG. The transmission source 102 generates a sound wave or a vibration wave. The generated sound waves or vibration waves are transmitted to the measurement target area 101. The receiver 103 measures a sound wave or a vibration wave generated from the transmission source 102 at each point of the measurement target area 101. The characteristics of the sound waves or vibration waves measured by the respective receivers 103 may vary depending on the underground structure of the area 101 to be measured.

The signal processing device 104 processes the sound wave or vibration wave measured by the receiver 103 to generate image data for the measurement target area 101. [ For generating image data, the signal processing device 104 may generate a velocity model for the region 101 to be measured. The velocity model can represent the velocity distribution of the acoustic wave with respect to the measurement target area 101. [ For example, the acoustic wave velocity in the measurement target area 101 varies depending on the characteristics of the respective points inside the measurement target area 101. The characteristics of each point are, for example, the type or density of the medium. When the velocity model representing the velocity distribution of the elastic wave is obtained, the underground structure of the measurement target region 101 can be imaged from the velocity model.

The signal processing device 104 is capable of generating a velocity model through waveform inversion. Waveform inversion is to create an initial velocity model for the region of interest and to update the initial velocity model repeatedly using measurements obtained in the same region of interest. That is, the initial velocity model is set as the reference velocity model to update the reference velocity model. Through the above-described repetitive updating of the speed model, a speed model having characteristics similar to the actual characteristics of the area of interest can be obtained.

For example, the signal processing device 104 sets a velocity model of the measurement target area 101, obtains seismic data from the measurement target area 101, and then sets the velocity model of the measurement target area 101 using the acquired elastic wave data A speed model similar to the velocity distribution of the actual seismic waves of the measurement target area 101 is obtained by repeatedly updating the velocity model.

2 shows the principle of waveform inversion according to an embodiment of the present invention.

Referring to Fig. 2, the measurement target region 201 has a certain characteristic V. Fig. The characteristic V may be various characteristics such as the velocity, density, or temperature of the elastic wave with respect to each part of the measurement target area 201. If an input S is applied to the region 201 to be measured, the output d can be observed according to the characteristic V thereof. At this time, the output d depends on the input S and the characteristic V. For example, even if the same input S is given, if the characteristic V of the measurement target region 201 is changed, the output d thereof can also be changed.

In Fig. 2, assuming that the characteristic V is the velocity characteristic of the region 201 to be measured, the region 201 to be measured can be modeled as a parameter m regarding the velocity. The parameter m regarding the velocity may be set to have a uniform velocity distribution at the beginning, since it is updated repeatedly thereafter. When the speed parameter m is set, the output u when the virtual input S is applied to the modeled measurement area 202 can be obtained. That is, the output u corresponds to the output d, and if the velocity parameter m is set equal to the actual characteristic V, the output u will have the same value as the output d. It is possible to obtain the same parameter m as the actual characteristic V by adjusting the speed parameter m so that the output u becomes equal to the output d.

Through such waveform inversion, the elastic-wave imaging apparatus using the repetitive application of the direct waveform inversion according to the present embodiment obtains the measurement data d and u as the modeling data, and obtains the velocity data such that the difference between the measurement data d and the modeling data u is minimized After estimating the characteristic V for the measurement target area by adjusting the parameter m, the velocity model and the image data can be generated using the obtained velocity parameter m or the estimated characteristic V. [

FIG. 3 shows a detailed configuration of an elastic-wave imaging apparatus using repetitive application of direct waveform inversion according to an embodiment of the present invention.

In one aspect, the underground structure imaging device includes a source, a source, and a signal processing device. A detailed description of the transmission source, the reception source, and the signal processing apparatus has been described above, and the signal processing apparatus shown in FIG. 1 includes the waveform inversion unit 300 and the underground structure imaging unit 310.

An elastic wave imaging apparatus using repetitive application of direct waveform inversion imaging an underground structure of a measurement target area according to one aspect includes a waveform inversion unit 300 and an underground structure imaging unit 310.

In one embodiment, the underground structure imaging unit 310 images the underground structure from the updated reference speed model. The underground structure imaging unit 310 outputs a signal received from the waveform inversion unit 300. The underground structure imaging unit 310 is, for example, a display device for outputting data processed in an image or graph format.

In one embodiment, the waveform inversion unit 300 updates the reference velocity model while moving the frequency band in a predetermined order using a virtual scattering source and a parameter perturbation calculated from the updated reference wave field, And updates the reference speed model updated in the previous frequency band to the reference speed model in the next frequency band.

The process of updating the reference velocity model using a virtual scattering source and a parameter perturbation calculated from the updated reference wave field will be described later.

The waveform inversion unit 300 may create an initial velocity model for the region of interest and iteratively update the initial velocity model using the measurements obtained in the same region of interest. That is, the initial velocity model and the reference velocity model may be the same in order to update the reference velocity model by setting the initial velocity model as the reference velocity model.

In one embodiment, the waveform inversion unit 300 sets the updated reference speed model in the last frequency band back to the reference speed model in the first frequency band, and updates the updated speed model in the previous frequency band to the next frequency band The reference speed model is set and updated. After updating the speed model for the first frequency band and moving to the last frequency band, the updated speed model is again assumed to be the reference speed model of the first frequency band and the inverse calculation is performed to obtain the final speed model close to the actual medium. For example, when updating the speed model in ascending order of frequency, the waveform inversion section 300 updates the speed model in the largest frequency band, which is the last frequency band, and updates the updated speed model to the reference speed model of the lowest frequency band To update the speed model. After updating the velocity model for the lowest frequency band, the frequency strategy of updating the velocity model can be performed more than once.

The above-described update of the speed model is not limited to the update in ascending order of frequency, but the speed model can be updated in descending order or in any order.

In one embodiment, when the frequency is a complex value, the waveform inversion unit shifts the frequency band while changing the damping constant for adjusting the imaginary part and the size of the real part in the set order. If we denote the frequency f of Hertz unit as angular frequency w (omega) with unit rad / s, we have w = 2 * pi * f without damping, where pi is 3.14159, that is, the ratio of the circumference. For complex angular frequencies including damping, w = 2 * pi * f + i * alpha. i is the square root of the unit complex number (-1) sqrt (-1) and alpha is the damping constant representing the damping strength. That is, the frequency is a real number, an imaginary number, or a complex number, and the waveform inversion section updates the velocity model by shifting the frequency band in the form of a real number, an imaginary number, and a complex number. If the frequency is a complex number, the waveform inversion unit updates the velocity model by shifting the frequency band while changing the damping constant represented by alpha and the magnitude f of the real part in the order set respectively.

In one embodiment, the set order is in ascending order, descending order, or any order. When the damping constant is 0, that is, when the frequency is a real number, the speed model is updated while changing only the size of the real part according to the set order.

When updating the velocity model in ascending order of real part size, the waveform inversion unit 300 uses a sequential frequency strategy to update the velocity model by moving to the higher frequency band after updating the velocity model for the lower frequency band. The waveform inversion unit 300 updates the initial velocity model by updating the velocity model for the lowest frequency band and setting the updated velocity model as the initial velocity model for the second lowest frequency band. The initial velocity model is the reference velocity model.

For example, if the selected frequency range is, for example, 15 Hz at 4 Hz and the selected frequency range is 0.2 Hz, the speed model is updated in the order of 4.0 Hz, 4.2 Hz, 4.4 Hz, ..., 15 Hz.

The method of updating the speed model is not limited to the ascending order described above and the speed model can be updated by setting the speed model updated in the previous frequency band in descending order of frequencies or in an arbitrary order to the initial speed model of the next frequency band . This iterative calculation is performed until the speed model is updated in the last frequency band.

The case where the damping constant is not 0 is as follows.

(4.0Hz, 0.2), (4.0Hz, 0.4) when the damping constants change from 0.2 to 1.0 in increments of 0.2 and the real part changes from 0.2 to 4 Hz in the ascending order from 4 Hz to 15 Hz, (4.0 Hz, 0.6), 4.0 Hz, 0.8, 4.0 Hz, 1.0, 4.2 Hz, 0.2, 4.2 Hz, 0.4, 4.2 Hz, 0.6, (15 Hz, 1.0), (15 Hz, 0.2), (15 Hz, 0.4), (15 Hz, 0.6) The velocity model is updated in the band.

An example in which the damping constants and the real part are in the same range and the same range as in the above example but the damping constant is descending and the real part is in ascending order is as follows.

(4.0 Hz, 1.0), (4.0 Hz, 0.8), (4.0 Hz, 0.6), (4.0 Hz, 0.4), (15 Hz, 0.6), (15 Hz, 0.4), (15 Hz, 0.4), (4.2 Hz, , 0.2), and the velocity model is updated in each frequency band.

An example in which the damping constants and the real part are in the same range and the same range as in the above example but the damping constant is descending and the real part is descending is as follows.

(15 Hz, 1.0), (15 Hz, 0.8), (15 Hz, 0.6), (15 Hz, 0.4) (4.0Hz, 0.4), (4.2Hz, 0.2), (4.0Hz, 1.0), (4.0Hz, 0.8), (4.0Hz, 0.6) , 0.2), and the velocity model is updated in each frequency band.

The damping constants and the real part are in the same range and the same range as in the above example, but the damping constants are in the descending order of two times and the real part is in ascending order.

(15Hz, 1.0), (15Hz, 1.0), (15Hz, 0.8), (15Hz, 0.8), ..., (4.2Hz, 0.4), (4.2Hz, (4.0 Hz, 0.6), (4.0 Hz, 0.6), 4.0 (4.0 Hz, 1.0), 4.0 Hz, 0.8, 4.0 Hz, Hz, 0.4), (4.0 Hz, 0.4), (4.0 Hz, 0.2), (4.0 Hz, 0.2), and the velocity model is updated in each frequency band.

The method of shifting the frequency band is not limited to the above example, and the real part size and the damping constant of the frequency can be set to various intervals and various ranges. The above-described flowcharts are not limited to the descending order and the ascending order, but may be in any order.

In one aspect, the waveform inversion section 300 includes a reference wave length calculation section 301, a virtual scatter transmission source calculation section 302, a scattered wave length calculation section 303, a reference wave length update section 304, A perturbation calculation unit 305 and a reference velocity model updating unit 306. [

In one embodiment, the reference wave length calculating unit 301 calculates a reference wave length, which is a solution of the wave equation, using data collected from at least one of the transmitter and the receiver. The collected data includes, for example, p wave velocity, wave field, source function, and the like.

In one embodiment, the reference wave length calculating unit 301 calculates a reference wave length by using a finite element method or a finite difference method, using the Fourier transform to convert the frequency equation of the frequency domain converted from the waveform equation of the time domain, . The process of calculating the reference wave field by the reference wave field calculation unit 301 is as follows.

The acoustic approximation of the 2D wave equation in the time domain can be given as Equation (1).

(One)

In Equation (1)

Is the p-wave velocity,

Is a pressure wave field,

Is a source wavelet function in the time domain.

The Fourier transform of equation (1) yields equation (2) which is a wave equation in the frequency domain.

(2)

In Equation (2)

Is the temporal angular frequency,

Is a pressure wave field in the frequency domain,

Is a source wavelet function in the frequency domain. Forward modeling of wave propagation simulation is a numerical solution of a wave form equation using a spatial approximation such as a finite difference method or a finite element method. Based.

Equation (2) can be rewritten as Equation (3) using the finite-difference method in the spatially domain.

(3)

Wow

Represent velocity and frequency domain wavefields at the grid point x, z of the velocity model, respectively.

Represents the lattice spacing in the model. In the frequency domain, the numerical solution of the waveform equation is obtained by solving the equation in matrix form. It is known that equation (2), which is a wave equation in the frequency domain, can be expressed by a general determinant as in equation (4).

(4)

S is the complex impedance matrix of the velocity model, u is the wave field and f is the source vector in the frequency domain.

(5)

(6)

(7)

The complex impedance matrix S is an mxm band matrix, and m is the number of lattice points in the velocity model. Assuming that the matrix S and S ₀ are the complex impedance matrix of the actual velocity model for the waveform inversion and the complex impedance matrix of the reference velocity model,

Can be represented by two matrices.

(8)

(9)

Is an impedance difference matrix

Lt; / RTI >

Is the P wave velocity of the Kth lattice point in the real velocity model

Is the P wave velocity of the Kth lattice point in the reference velocity model. Diagonal

Is the parameter between the velocity of the actual medium and the velocity of the reference medium (=

) Is proportional to the difference.

The main purpose of direct waveform inversion is to obtain the parameter perturbation directly between the real velocity model and the reference velocity model. The wave field (U ₀ ) of the reference velocity model is obtained by solving the following waveform equation (10) in matrix form.

(10)

U ₀ is the reference wave field vector of the reference velocity model when the source vector of equation (10) is _equal to the source vector of equation (4).

In one embodiment, the scattered wave field calculator 303 calculates the scattered wave field from the virtual scattered transmission source. The scattered wave field is the difference between the actual wave field and the reference wave field

Is expressed.

The equation (11) can be calculated from the difference between the equations (4) and (10).

(11)

Is an unknown matrix obtained through direct waveform inversion. Since the wave field is obtained at the receiver from field exploration,

It is an unknown vector. It is necessary to introduce the concept of a virtual scattering source to obtain a scattered wavefield based on a single scattering assumption and a Lippmann Shwinger equation. In Equation (11), the virtual scattering transmission source is described as Equation (12).

(12)

In the reference velocity model, the propagation of the virtual scattering source is the scattering of the scattered wave field to the wave field difference between the actual wave field and the reference wave field

).

In one embodiment, the virtual scattering source calculating unit calculates a virtual scattering source from the residual data, which is the difference between the actual wave field and the reference wave field. In one embodiment, the virtual scattering transmission source calculation unit 302 calculates the virtual scattering transmission source using the residual data and the function of green in the formula derived from the Lippmann Schwinger equation.

The virtual scattering transmission source is obtained from the data residual of Eq. (13).

(13)

Wow

Denote the receiver node and the total lattice point, respectively.

Wow

Represents observed field data and modeled data at the receiver node. Therefore, the left side of equation (13) is the same as the residual data.

Represents the Green's function at the receiver node (x ') when the impulsive source is located at an arbitrary point at the entire lattice point of the velocity model.

Virtual scattering source

Is a vector that emits a scattered wave field at the entire lattice point including the receiver node. Equation (13) is derived from the Lefman-Schwinger equation of scattering theory with a single scattering assumption.

The Green's function in Eq. (13) can be computed numerically from the reference velocity model with a single scattering assumption. As a numerical solution of Green's function, the inverse of the complex impedance of the reference velocity model (

) Can be used. Equation (13) can be expressed as Equation (14) using Equation (12).

(14)

The total number of grid points is m and the number of receiver nodes is n. In Equation (14)

And the other equation (12)

Is the wave field difference between the actual wave field and the reference wave field at the entire lattice point. The field data observed in the receiver position is only available in the waveform inversion process of the actual field data set. In Equation (14)

Is an mxm projection matrix that extracts a value corresponding to the receiver node at the entire lattice point. The projection matrix can be expressed by Equation (15).

(15)

Represents an nxn identity matrix. The projection matrix depends on the grid point indexing and assumes that the nth grid point from the first corresponds to the surface receiver node. The right side of equation (14)

Is a residual data vector as a difference between the field data and the model data. A least-squares method is introduced to derive a virtual scattering source vector from equation (14) and derive equation (16) as a normal equation.

(16)

The superscript * denotes the conjugate transpose of the matrix.

Is a symmetric matrix, i. E.

to be. Equation (16) can be rewritten as Equation (17).

(17)

- denotes the complex conjugate of the matrix. Eq. (17) is solved and the virtual scattering source vector at each frequency and each shot-gather

Can be obtained. Since the whole matrix is needed to calculate the virtual scattering source through Eq. (17), we can use a multifrontal matrix solver and a generalized minimal residual (Gmres) routine for numerical implementation Can be used.

In one embodiment, the reference wave field updating unit 304 updates the reference wave field by adding the reference wave field and the scattered wave field.

In the reference velocity model, the propagation of the virtual scattering source vector generates a scattered wavefield at the whole lattice point. The total wave field is the reference wave field

And scattered wave field

Can be added. The total wave field is the updated reference wave field. Because of the single scattering assumption, the generated wave field is not the same as the actual wave field propagated in the real velocity model, but the actual wave field can be replaced with the generated wave field. Because inverse waveform inversion utilizes parameter perturbation and single scattering assumption, the inverse result depends on the initial setting of the velocity model.

In one embodiment, the parameter perturbation calculation unit 305 calculates a parameter perturbation from the virtual scattering transmission source and the updated reference wave field. In one embodiment, the parameter perturbation calculator 305 calculates a parameter perturbation by applying a Newton method or a least squares method to the objective function.

Parameter perturbation

The parameter perturbation is obtained from the virtual scattering source and the generated wave field. The impedance difference matrix including the parameter perturbation in the diagonal elements between the impedance matrix of the actual velocity model and the impedance matrix of the reference velocity model can be calculated.

Is a diagonal factor proportional to the parameter perturbation between the real velocity model and the reference velocity model. Equation (11) is transformed into Equation (18) and an impedance difference matrix

Can be calculated.

(18)

The sum of the reference wave field and the scattered wave field

Is a recalculated wave field. The actual wave field at the whole lattice point can be replaced by the re-calculated wave field at the whole lattice point. This substitution allows direct inversion to overcome the limitations of seismic data stored solely in the receiver position from field exploration.

The direct waveform inversion method uses the

To provide parameter perturbations and to update the reference velocity model. In Equation (18), the impedance difference matrix

Virtual scattering source

And the re-

/ RTI > Complex number

(18) is a block system which is a component solved by a component since the components of the block are independent of each other.

Newton's method or least squares method can be used to minimize the residual between the left and right sides of equation (18). Virtual scattering source

And regenerated wave field

Is obtained in each of the shot gatherings, the objective function is the sum of the residuals of the various shot generators.

The objective function for the calculation is given by Eq. (19).

(19)

Residual

Is expressed by equation (20).

(20)

,

And

Lt; RTI ID = 0.0 > k < / RTI >

,

&Lt; / RTI > When the direct waveform inversion method is performed in the frequency domain, the variable may be a complex number.

(21)

(22)

(23)

The impedance difference matrix

It can be expressed as a complex number. The parameter perturbation may be a real number to correspond to the difference between the actual speed and the reference speed.

The real part of the reference velocity model may be needed to update the reference velocity model. Newton method or least squares method can be applied to Eq. (19) to obtain Eq. (24) which updates the parameter perturbation at the kth lattice point of the velocity model.

(24)

The Hessian matrix is shown in Eq. (25).

(25)

The gradient vector is the same as equation (26).

(26)

The equation (24) can be converted into the equation (27).

(27)

On the right side of equation (27)

and

Is converged to zero, so the resultant equation consists solely of c, d, e, and g.

The

Wow

Can be expressed as a complex number using.

(28)

Impedance difference matrix

The reference velocity model can be updated after the diagonal element of equation (28) is obtained from equation (28). Because the parameter perturbation is a mistake,

Is used to update the reference speed.

In one embodiment, the reference velocity model update unit 306 may update the reference velocity model from the data collected from at least one of the transmitter and the receiver, and the calculated parameter perturbation.

Equation (29) shows how to update the reference velocity model.

(29)

The direct waveform inversion can generate the virtual scattering source from the residual data by calculating the model data in the reference velocity model through Equation (10) and Equation (17). Direct waveform inversion can calculate the scattered wave field by propagating the virtual scattered transmission source in the reference velocity model, and reproduce the entire wave field by the sum of the reference wave field and the scattered wave field. The direct waveform inversion can obtain the parameter perturbation from the virtual scattered transmission source and the recalculated wave field via equation (28) and update the reference velocity model from equation (29). The equation (29) is obtained by updating the velocity model by using a sloth as a parameter, but is not limited thereto. The parameters can be various parameters such as the square of the velocity, the reciprocal of the velocity, and so on.

4 is a flowchart of an elastic wave imaging method using repetitive application of direct waveform inversion according to an embodiment.

In order to image the underground structure of the measurement target area according to one aspect, the elastic wave imaging method using repetitive application of direct waveform inversion may include a waveform inversion step (S410) and an underground structure imaging step (S420).

In one embodiment, the underground structure imaging (S420) images the underground structure from the updated reference velocity model. In the underground structure imaging step S420, the signal received from the waveform inversion unit 300 is output. In the underground structure imaging step S420, the display device can output data processed in an image or a graph format.

In one embodiment, the waveform inversion step S410 is performed by moving a frequency band in a predetermined order using a virtual scattering source and a parameter perturbation calculated from the updated reference wave field, (S411 to S416) of updating the reference speed model updated in the previous frequency band to the reference speed model of the next frequency band will be described later.

The waveform inversion step S410 may create an initial velocity model for the region of interest and may iteratively update the initial velocity model using measurements obtained in the same region of interest. That is, the initial velocity model and the reference velocity model may be the same in order to update the reference velocity model by setting the initial velocity model as the reference velocity model.

In one embodiment, the waveform inversion step S410 sets the updated reference speed model in the last frequency band back to the reference speed model in the first frequency band, and updates the updated speed model in the previous frequency band to the next frequency band The reference speed model is set and updated. After updating the speed model for the first frequency band and moving to the last frequency band, the updated speed model is again assumed to be the reference speed model of the first frequency band (S417, S418, S419) Can be obtained. For example, when updating the speed model in ascending order of frequency, the waveform inversion step S410 updates the speed model in the largest frequency band, which is the last frequency band, and updates the updated speed model to the reference speed model of the lowest frequency band To update the speed model (S418, S419). After updating the velocity model for the lowest frequency band, the frequency strategy of updating the velocity model can be performed more than once.

In one embodiment, when the frequency is a complex value, the waveform inversion step S410 shifts the frequency band while changing the damping constant for adjusting the imaginary part and the size of the real part in the set order. If we denote the frequency f of Hertz unit as angular frequency w (omega) with unit rad / s, we have w = 2 * pi * f without damping, where pi is 3.14159, that is, the ratio of the circumference. For complex angular frequencies including damping, w = 2 * pi * f + i * alpha. i is the square root of the unit complex number (-1) sqrt (-1) and alpha is the damping constant representing the damping strength. That is, the frequency is a real number, an imaginary number, or a complex number, and the waveform inversion section updates the velocity model by shifting the frequency band in the form of a real number, an imaginary number, and a complex number. If the frequency is a complex number, the waveform inversion unit updates the velocity model by shifting the frequency band while changing the damping constant represented by alpha and the magnitude f of the real part in the order set respectively.

When the speed model is updated in ascending order of the real part, the waveform inversion step (S410) uses a sequential frequency strategy to update the speed model while moving to the higher frequency band after updating the speed model for the lower frequency band. The waveform inversion step S410 updates the initial velocity model by updating the velocity model for the lowest frequency band and setting the updated velocity model as the initial velocity model for the second lowest frequency band.

The case where the damping constant is not 0 is as follows.

5 is a detailed flowchart of the waveform inversion step S410 according to an embodiment.

In one aspect, the waveform inversion step S410 includes a reference wave length calculation step S411, a virtual scatter transmission source calculation step S412, a scattered wave length calculation step S413, a reference wave length update step S414, A perturbation calculation step S415 and a reference velocity model update step S416.

In one embodiment, the reference wave length calculation step S411 calculates a reference wave length, which is a solution of the wave equation, using data collected from at least one of the transmitter and the receiver. The collected data includes, for example, p wave velocity, wave field, source function, and the like.

In one embodiment, the reference wave length calculation step S411 is a step of calculating a reference wave length field by using a finite element method or a finite difference method, using a Fourier transform to convert a frequency equation of a frequency domain from a time domain waveform equation, . The reference wave length calculation step S411 calculates the reference wave length as follows.

(One)

In Equation (1)

Is the p-wave velocity,

Is a pressure wave field,

Is a source wavelet function in the time domain.

(2)

In Equation (2)

Is the temporal angular frequency,

Is a pressure wave field in the frequency domain,

(3)

Wow

(4)

(5)

(6)

(7)

Can be represented by two matrices.

(8)

(9)

Is an impedance difference matrix

Lt; / RTI >

Is the P wave velocity of the Kth lattice point in the real velocity model

) Is proportional to the difference.

(10)

In one embodiment, the scattered wave field computing step S413 produces a scattered wave field from the virtual scatter transmission source. The scattered wave field is the difference between the actual wave field and the reference wave field

Is expressed.

(11)

It is an unknown vector. We introduce the concept of a virtual scattering source to obtain a scattered wavefield based on the single scattering assumption and the Lippmann-Shwinger equation. In Equation (11), the virtual scattering transmission source is described as Equation (12).

(12)

).

In one embodiment, the virtual scattering source calculating step calculates a virtual scattering source from the residual data, which is the difference between the actual wave field and the reference wave field. In one embodiment, the virtual scatter transmission source calculation step S412 calculates a virtual scatter transmission source using a function of residual data and Green in an expression derived from a Lippmann Schwinger equation.

(13)

Wow

Denote the receiver node and the total lattice point, respectively.

Wow

Virtual scattering source

(14)

And the other equation (12)

(15)

(16)

The superscript * denotes the conjugate transpose of the matrix.

Is a symmetric matrix, i. E.

to be. Equation (16) can be rewritten as Equation (17).

(17)

^- denotes the complex conjugate of the matrix. Eq. (17) is solved and the virtual scattering source vector at each frequency and each shot-gather

In one embodiment, the reference wave field updating step S414 updates the reference wave field by adding the reference wave field and the scattered wave field.

And scattered wave field

Can be added. The total wave field is the updated reference wave field. The actual wave field can be replaced with the generated total wave field. Because inverse waveform inversion utilizes parameter perturbation and single scattering assumption, the inverse result depends on the initial setting of the velocity model.

In one embodiment, the parameter perturbation calculation step S415 calculates a parameter perturbation from the virtual scatter transmission source and the updated reference wave field. In one embodiment, the parameter perturbation calculation step S415 calculates the parameter perturbation by applying a Newton method or a least squares method to the objective function.

Parameter perturbation

Can be calculated.

(18)

The sum of the reference wave field and the scattered wave field

The direct waveform inversion method uses the

Virtual scattering source

And the re-

/ RTI > Complex number

The Newton method or the least squares method can be used to minimize the residual between the left and right sides of equation (18). Virtual scattering source

And regenerated wave field

The objective function for the calculation is given by Eq. (19).

(19)

Residual

Is expressed by equation (20).

(20)

,

And

Lt; RTI ID = 0.0 > k < / RTI >

,

(21)

(22)

(23)

The impedance difference matrix

(24)

The Hessian matrix is shown in Eq. (25).

(25)

The gradient vector is the same as equation (26).

(26)

The equation (24) can be converted into the equation (27).

(27)

On the right side of equation (27)

and

The

Wow

Can be expressed as a complex number using.

(28)

Impedance difference matrix

Is used to update the reference speed.

In one embodiment, the reference velocity model update step (S416) may update the reference velocity model from the data collected from at least one of the transmitter and the receiver and the computed parameter perturbations.

Equation (29) shows how to update the reference velocity model.

(29)

The direct waveform inversion can generate the virtual scattering source from the residual data by calculating the model data in the reference velocity model through Equation (10) and Equation (17). Direct waveform inversion can calculate the scattered wave field by propagating the virtual scattered transmission source in the reference velocity model, and reproduce the entire wave field by the sum of the reference wave field and the scattered wave field. The direct waveform inversion can obtain the parameter perturbation from the virtual scattered transmission source and the recalculated wave field via equation (28) and update the reference velocity model from equation (29).

Figure 6 shows the Marmousi velocity model.

In detail, FIG. 6 shows a Marmousi velocity model applying the direct waveform inversion method to verify the proposed direct waveform inversion method.

7 (a) shows a shot-gather seismogram and (b) shows a frequency spectrum of a Marmousi data set.

Referring to FIG. 7A, an arbitrary field data set may be generated by 192 shots of 48m intervals. Each shot gather includes 384 receivers at 24m intervals. A finite-difference scheme can be used, and an absorptive boundary condition is applied to the left, right and bottom of the model. Free surface conditions apply to the top node. Direct waveform inversion can be used in the frequency domain.

Referring to FIG. 7 (b), the arbitrary field data has a frequency spectrum ranging from 0 to 40 Hz. In the inversion process, 161 frequencies in the range of 4 to 20 Hz and 0.1 Hz intervals can be used, and the source wavelet signature can be estimated by the Newton method or the least square method in the frequency domain.

8 shows an initial velocity model.

Referring to FIG. 8, the Smoothed velocity structure of the Marmousi model is used as the initial velocity model.

Figures 9 (a), 9 (b) and 9 (c) each show a reference wavelet propagated at a shot position of 4.8 Km in an initial velocity model of 5.0 Hz, 10.0 Hz and 15.0 Hz.

9 shows the frequency domain wave field propagated in the initial velocity model, which is the reference wave field in equation (10). The inversion method calculates the velocity model perturbations directly from the virtual scattering source and the regenerated wave field. From the residual data using the Lehman-Schwinger equation, the virtual scattering source vector

Can be calculated.

Each of Figs. 10A, 10B and 10C shows a virtual scattering source obtained from residual data corresponding to a point source of a shot position of 4.8 Km at 5.0 Hz, 10.0 Hz and 15.0 Hz do.

Figure 10 shows the scattered wave field

0.0 > scattering < / RTI > A multifrontal matrix solver can be used to solve the forward modeling and GMRES routines to solve the virtual scattering source in Eq. (17). Because multiple front solvers can decompose the inverse system matrix and store it in RAM, multiple front solvers can solve multiple numbers of the right-side vector until the matrix decomposition is done.

Since multiplying the inverse matrix is equivalent to solving the system matrix, it is not necessary to store the inverse matrix of complex impedance to solve the virtual scattering source. The direct inversion method does not require much memory compared to the propagation inversion.

The direct inversion method requires less model update because the speed model is updated once per frequency. However, the calculation of the virtual scattering source requires a relatively long computation time with multiple front solvers and GRMES routines.

The computation time of the matrix decomposition by the multiple front solver in the inversion test is about 0.32 seconds and the computation time of the decomposed matrix for the right side vector by the system matrix is about 0.04 seconds. The time to solve the virtual scattering transmission source is about 11.96 seconds because of the system matrix that repeats Equation (17) about 200 times.

Figures 11 (a), (b) and (c) each show scattered wave filed generated by the virtual scattering source shown in Figure 10 at 5.0 Hz, 10.0 Hz and 15.0 Hz.

Referring to FIG. 11 and equation (12), the scattered wave field

Is obtained by propagating the virtual scattering source in the reference velocity model.

Each of Figs. 12 (a), (b) and (c) shows the generated wave at the shot position of 4.8 Km, which is the sum of the reference wave field of Fig. 9 and the scattered wave field of Fig. 10 at 5.0 Hz, 10.0 Hz and 15.0 Hz Field.

Referring to Fig. 12, the wave field propagated in the structure immediately below the surface is replaced by the recalculated wave field, which is the sum of the reference wave field and the scattered wave field.

Figures 13 (a), 13 (b) and 13 (c) each show the actual wave field at a shot position of 4.8 Km at 5.0 Hz, 10.0 Hz and 15.0 Hz.

Referring to Fig. 13, the re-calculated wave field is similar to the actual wave field.

Figures 14 (a), 14 (b) and 14 (c) show parameter perturbations at 5.0 Hz, 10.0 Hz and 15.0 Hz, respectively.

Referring to FIG. 14 and (28), the parameter perturbation is computed as a model perturbation using the virtual scatter transmission source and the recalculated wave field. The initial velocity model can be updated by adding the velocity model perturbation using equation (29). The process from source signature estimation to speed model update is performed at a single frequency. Before moving to the next frequency, the speed update is performed and the frequency is shifted continuously from low frequency to high frequency.

Figure 15 shows the inversed velocity model.

Referring to FIG. 15, an inverted velocity model is similar to an actual velocity model in stratified layers and structures and shows an anticline structure at the bottom center of an actual model.

16 (a), 16 (b), and 16 (c) are diagrams each showing the Shotgarder seismic data of synthesized field data, the synthesized seismic data in the initial velocity model, and the seismic data obtained by direct waveform inversion / RTI >

The resultant velocity model of the direct waveform inversion can be verified by any elastic wave recording generated in the inversion process. Time domain seismic recording is generated by inverse Fourier transform of the re-computed wave field in the frequency domain. Any acoustic wave record is the same as the propagation of the source estimated in the inversed velocity model, and any acoustic wave record matches the corresponding shotgarder acoustic wave record of the field data.

17 (a) and 17 (b) show Shotgarder seismic data of a noise-mixed Marmousi data set with signal-to-noise ratios of 12.75 and 25.50, respectively.

To check the sensitivity of the inverse result in incoherent noise, you can test the direct waveform inversion of the noise-containing data. Referring to Figures 2 and 12, Gaussian random noise is added to the Marmousi data set with a signal to noise ratio of 12.75 and 25.50. Due to strong incoherent noise, noise data (noisy data) with a signal-to-noise ratio of 12.75 can not account for significant reflections or refraction signals. Direct waveform inversion is performed under the same conditions as the Marmousi data set test.

18 (a) and 18 (b) show the inverse model of the noise in the noise-mixed Marmousi data set with S / N ratios of 12.75 and 25.50, respectively.

Figure 18 shows the inverse result. Since the important signals of the data include random noise, the waveform inversion result is affected by the amplitude of the noise data. FIG. 13 shows that the inverse result of the data having the signal-to-noise ratio of 25.50 is stratified to the anticline structure of the data having the signal-to-noise ratio of 12.75.

Since the direct waveform inversion method uses residual data to generate a virtual scattered transmission source and a scattered wave field, the result of the direct waveform inversion depends on the quality of the input data such as the amplitude of the signal. If input seismic data contain noise, the virtual scattering source can be made to regenerate the scattered wave field containing the noise signal of the input data.

Since the direct waveform inversion method uses a solution of the propagation equation with forward modeling in the whole process like propagation inversion, the multiple in the input data need not be removed. Coherent noise must be removed and otherwise direct inversion is induced to produce a virtual scattering source for the scattered wave field containing coherent noise. If the input data has a sufficiently high signal-to-noise ratio to preserve coherent noise and significant signals, then direct waveform inversion can provide a reasonable solution to the velocity model.

Direct waveform inversion is a method based on perturbation theory that measures the p-wave velocity of a medium below the surface from seismic exploration data. The inversion method uses residual data in a manner similar to propagation inversion, but does not use the steepest descent method to reduce residual data. The computational expense required to perform the direct inversion is much smaller than the propagation inversion based on the most rapid method requiring many forward modeling and backpropagation steps. Since the direct waveform inversion directly inverses the parameter perturbation, it does not require a step length to update the speed model, but it can also adjust the speed update amount by multiplying it by a value smaller than 1.

The inversion process generates a model acoustic wave record similar to the acoustic wave record of the field data corresponding to the shot gathers. Any shot-girder seismic records generated in the inversion process demonstrate that the inverse velocity model provides a reasonable solution.

In real field data applications, noise must be removed from seismic exploration data. Most intermediate outputs are associated with a virtual scatter source from the residual data. The amplitude of the signal in the input data can have a significant effect on the inverse result. The test results of the noise-containing data show that the low signal-to-noise ratio of the input data leads to inaccurate estimation of the parameter perturbation. In order to obtain the proper inversion velocity model, the noise signal, especially the coherent noise, in the input signal must be removed without ruin of important signals such as reflections and refractions.

While the present invention has been particularly shown and described with reference to exemplary embodiments thereof, it is clearly understood that the same is by way of illustration and example only and is not to be taken by way of limitation and that those skilled in the art will recognize that various modifications and equivalent arrangements may be made therein. It will be possible. Accordingly, the true scope of protection of the present invention should be determined only by the appended claims.

101: Measurement area
102: source
103: Receiver
104: Signal processing device
201: Measuring area
202: Modeled measurement area
300: waveform inversion section
301: reference wave length calculating section
302: virtual scattering transmission source calculating unit
303: scattered wave field calculation section
304: reference wave length renewal unit
305: Parameter fluctuation calculating section
306: Reference speed model update unit
310: Underground structure visualizer

Claims

1. An elastic-wave imaging apparatus using repetitive application of direct waveform inversion imaging an underground structure of a measurement target area,
The reference velocity model is updated while shifting the frequency band in a predetermined order using a virtual scattering source and a parameter perturbation calculated from the updated reference wave field, To a reference velocity model of the next frequency band to update the reference velocity model; And
An underground structure imaging unit for imaging an underground structure from an updated reference velocity model; , &Lt; / RTI &
The waveform inversion section updates the reference velocity model in the last frequency band and sets the reference velocity model in the first frequency band again to update the reference velocity model while repeating the frequency band in the set order. The reference speed model of the next frequency band is set and updated
Elastic wave imaging device using repetitive application of direct waveform inversion.

delete

The apparatus of claim 1, wherein the waveform inversion section
When the frequency is a complex value, the damping constant for adjusting the imaginary part and the magnitude of the real part are changed in the set order and the frequency band is shifted
Elastic wave imaging device using repetitive application of direct waveform inversion.

4. The method of claim 1,
In ascending, descending, or arbitrary order.
Elastic wave imaging device using repetitive application of direct waveform inversion.

The apparatus of claim 1, wherein the waveform inversion section
A reference wave length calculation unit for calculating a reference wave length, which is a solution of the wave equation, using data collected from at least one of a transmitter and a receiver;
A virtual scattering transmission source calculating unit for calculating a virtual scattering transmission source from residual data that is a difference between an actual wave field and a reference wave field;
A scattered wave field calculation unit for calculating a scattered wave field from a virtual scattering source;
A reference wave length updating unit for updating the reference wave length by adding the reference wave length and the scattered wave length;
A parameter perturbation calculation unit for calculating a parameter perturbation from the virtual scattering transmission source and the updated reference wave field;
A reference velocity model updating unit for updating the reference velocity model from the data collected from at least one of the transmitter and the receiver and the calculated parameter perturbation; To
Elastic wave imaging device using iterative application of direct waveform inversion including.

The apparatus of claim 5, wherein the reference wave length calculation unit
A frequency domain waveform equation transformed from a time domain waveform equation is transformed into a reference wave field by a finite element method or a finite difference method using a Fourier transform
Elastic wave imaging device using repetitive application of direct waveform inversion.

6. The apparatus according to claim 5, wherein the virtual scattering transmission source calculating unit
In the formula derived from the Lippmann Schwinger equation, the residual data and the function of green are used to calculate the source of the virtual scattering source
Elastic wave imaging device using repetitive application of direct waveform inversion.

6. The apparatus of claim 5, wherein the parameter perturbation calculating unit
The Newton method or the Least Squaring method is applied to the objective function to calculate the parameter perturbation
Elastic wave imaging device using repetitive application of direct waveform inversion.

A method for imaging an elastic wave using repetitive application of direct waveform inversion imaging an underground structure of a measurement target area,
The reference velocity model is updated while shifting the frequency band in a predetermined order using a virtual scattering source and a parameter perturbation calculated from the updated reference wave field, To a reference velocity model of the next frequency band and updates the reference velocity model; And
An underground structure imaging step imaging the underground structure from the updated reference velocity model; , &Lt; / RTI &
In the waveform inversion step, the reference speed model updated in the last frequency band is set again as the reference speed model of the first frequency band, and the reference speed model is updated while the frequency band is shifted in the set order repeatedly. Is set to the reference speed model of the next frequency band and is updated
A Seismic Imaging Method Using Repeated Application of Direct Wave Inversion.

delete

10. The method of claim 9, wherein the waveform inversion step
When the frequency is a complex value, the damping constant for adjusting the imaginary part and the magnitude of the real part are changed in the set order and the frequency band is shifted
A Seismic Imaging Method Using Repeated Application of Direct Wave Inversion.

10. The method of claim 9,
In ascending, descending, or arbitrary order.
A Seismic Imaging Method Using Repeated Application of Direct Wave Inversion.

10. The method of claim 9, wherein the waveform inversion step
A reference wave length calculation step of calculating a reference wave length, which is a solution of the wave equation, using data (p wave velocity, wave field, transmission source function) collected from at least one of a transmitter and a receiver;
A virtual scattering transmission source calculating step of calculating a virtual scattering transmission source from residual data that is a difference between an actual wave field and a reference wave field;
A scattered wave field computing step of computing a scatter wave field from a virtual scattering transmission source;
A reference wave field updating step of updating a reference wave field by adding a reference wave field and a scattered wave field;
A parameter perturbation calculation step of calculating a parameter perturbation from the virtual scattering transmission source and the updated reference wave field;
A reference velocity model updating step of updating the reference velocity model from the data collected from at least one of the transmitter and the receiver and the calculated parameter perturbation; To
A method of elastic wave imaging using repetitive application of direct waveform inversion including.

14. The method of claim 13, wherein the reference wave length calculation step
A frequency domain waveform equation transformed from a time domain waveform equation is transformed into a reference wave field by a finite element method or a finite difference method using a Fourier transform
A Seismic Imaging Method Using Repeated Application of Direct Wave Inversion.

14. The method of claim 13, wherein the virtual scattering transmission source calculating step
In the formula derived from the Lippmann Schwinger equation, the residual data and the function of green are used to calculate the source of the virtual scattering source
A Seismic Imaging Method Using Repeated Application of Direct Wave Inversion.

14. The method of claim 13, wherein the parameter perturbation calculation step
A method of elastic wave imaging using iterative application of direct waveform inversion to calculate parameter perturbations by applying Newton method or least squares method to objective function.