The invention is described more fully hereinafter with reference to the accompanying drawings, in which exemplary embodiments of the invention are shown. This invention may, however, be embodied in many different forms and should not be construed as limited to the exemplary embodiments set forth herein. Rather, these exemplary embodiments are provided so that this disclosure is thorough, and will fully convey the scope of the invention to those skilled in the art. Like reference numerals in the drawings denote like elements.

Waveform inversion according to an exemplary embodiment of the present invention is a process which deduces information (for example, velocity model or density model of a target region that is to be inverted) about subsurface structure based on the measured data in a field. The waveform inversion can include a modeling process in which an analyzer sets up an initial velocity model and obtains a modeled value for the subsurface velocity model.

For example, when subsurface structure is updated using waveform inversion according to the exemplary embodiment, by comparing modeled values computed through modeling with measured values obtained through actual field exploration to create a new subsurface velocity model, and iteratively comparing modeled values of the new subsurface velocity model with the measured values obtained through actual field exploration to minimize a misfit there between, an subsurface velocity model which is close to the true subsurface velocity model can be obtained.

The waveform inversion according to the exemplary embodiment can be implemented by a computation apparatus for processing a variety of signals to generate image data for visualizing subsurface structure of a target region, a computer-readable recording medium on which signal processing algorithms are recorded, and a method for visualizing subsurface structure through the computation apparatus or the recording medium, etc.

A waveform inversion principle according to an exemplary embodiment of the present invention will be described with reference to FIG. 1, below.

In FIG. 1, V denotes physical characteristics of the subsurface structure, S denotes an input applied to V, and D denotes an output when S is applied to V.

The waveform inversion is aimed at finding the characteristic value V of the subsurface structure using a measured value D. If a velocity distribution (that is, a seismic wave velocity model) among the characteristics of the subsurface structure is obtained, the subsurface structure can be easily visualized using the velocity distribution. For this reason, it is assumed that the characteristic V is a velocity distribution, the input S is a source, and the output D is seismic data.

Measured seismic data Dreal can be actually obtained through actual field exploration of a target region, and estimated seismic data Dest can be obtained from estimated values (Vest, Sest) obtained by modeling of the target region. By updating initially estimated velocity data Vest and source data Sest and iteratively performing the updating process, a difference between the measured seismic data Dreal and the estimated seismic data Dest is minimized and an estimated velocity distribution Vest is equalized to an actual velocity distribution Vreal. Here, V representing subsurface structure can represent only a velocity distribution as described above or it can represent velocity/density, impedance, Lame constant/density distribution.

FIG. 2 is a schematic diagram for explaining a subsurface velocity structure visualization apparatus according to an exemplary embodiment of the present invention.

Referring to FIG. 2, the subsurface velocity structure visualization apparatus according to an exemplary embodiment of the present invention can include a source 102, a receiver 101, a signal processor 103, and a display unit 104.

The source 102 can generate a wave to a target region 105, and the receiver 101 can receive a seismic signal propagating through the target region 105. The signal processor 103 receives the seismic signal from the receiver 101, and processes the received seismic signal to generate image data so that the subsurface velocity structure of the target region 105 can be displayed through the display unit 104.

The signal processing performed by the signal processor 102 uses the waveform inversion described above, and the signal processor 103 performs the waveform inversion in the Laplace-Fourier domain. For example, the signal processor 103 performs domain transformation by multiplying a time-domain input by an exponential damping function e ^-st with a complex damping constant s and integrating the result of the multiplication with respect to time.

The Laplace-Fourier domain for the waveform inversion of the signal processor 103 will be described below in detail. The Laplace-Fourier transform for a time-domain wavefield can be expressed by Equation 1:

[Corrected under Rule 26 07.08.2009]
,(1)

where d(s) denotes data which is transformed to the Laplace-Fourier domain, d(t) denotes time-domain data, and s is a damping constant. The damping constant s is a complex number, which consists of a Laplace damping constantσ as a real part and an angular frequency ω as an imaginary part.

In Equation 1, if the damping constant s is real, the transformation is identical to the typical Laplace transform, and if the transformation is applied to the waveform inversion according to the exemplary embodiment, a zero-frequency component of a damped wavefield can be used.

In Equation 1, if the damping constant s is complex, the transformation is to a hybrid domain and is referred to as "the Laplace-Fourier transform" in this specification in order to distinguish it from the typical Laplace transform. That is, the Laplace-Fourier transform can be considered as an extension of the Laplace transform, and the waveform inversion in the Laplace-Fourier domain uses both low-frequency components and a zero-frequency component of a damped wavefield.

In order to compare a wavefield of the Laplace domain with a wavefield of the Laplace-Fourier domain, if it is assumed that a wavefield in the time domain is v(t), v(t) can be expressed by Equation 2:

[Corrected under Rule 26 07.08.2009]
(2)

If the time-domain wavefield v(t) is transformed using Equation 1, the transformed wavefield can be expressed by Equation 3:

[Corrected under Rule 26 07.08.2009]
(3)

It can be understood in Equation 3 that the Laplace-Fourier transformed wavefield v(s) is a function of an amplitude (a1 and b1) and phase (ω _c and ω _s ) of the time-domain wavefield v(t). For convenience of description, if it is assumed that b1 and ω _s are zero and ω _c is 20.0 or 40.0, the Laplace-Fourier transform v(t) can be expressed as shown in FIG. 3.

FIG. 4 shows a result obtained by the Laplace-transformation of the time-domain wavefield of FIG. 3, wherein s is a real number. As shown in FIG. 4, the amplitude of the wavefield increases as the Laplace damping constant increases, and the slopes of two curves are different according to the angular frequency of the wavefield.

FIG. 5 shows results obtained by Laplace-transforming the wavefield of FIG. 3, wherein s is a complex number. As shown in FIG. 5, unlike FIG. 4, both amplitude and phase are shown in a contour-like form.

As shown in FIG. 5, both the amplitude and the phase of the damped wavefield have peak spectrums of v(t) at 20 and 40, but also spread out to the nearby main frequency band. This indicates that when we use the damped wavefield of field data, which do not generate low-frequency components below 5 Hz for the waveform inversion, we can use "mirage-like" low-frequency components of the damped wavefield. Therefore, the waveform inversion in the Laplace-Fourier domain can obtain a velocity model including long-wavelength information using the low-frequency components.

FIG. 6 shows the Laplace-Fourier-domain wavefields when different Laplace damping constants are used. (a), (b) and (c) of FIG. 6 represent time-domain wavefields, and (d), (e) and (f) of FIG. 6 show frequency spectrums of Laplace-Fourier-domain wavefields, respectively. As shown in FIG. 6, as the Laplace damping constant increases, the amplitude and phase spectrums are smoothly attenuated. Particularly, the phase spectrum of 0 Hz to 5 Hz which is important in the waveform inversion notably changed with the Laplace damping constant.

Generally, waveform inversion has been performed in the time or frequency domain. However, since low-frequency components of field data are not reliable, there is a limitation in obtaining long wavelength information from field data. However, if waveform inversion is performed in the Laplace domain or the Laplace-Fourier domain as described above, low-frequency components near zero frequency can be used so that accurate information on subsurface velocity structure can be obtained.

If waveform inversion is performed in the Laplace domain, there is an advantage of obtaining long wavelength information using a zero-frequency component of the damped wavefield. However, since the penetration depth of the inversion is very sensitive to a maximum offset of seismic exploration data and a Laplace damping constant, the waveform inversion in the Laplace domain merely provides a velocity model in a smooth form.

Accordingly, in the current exemplary embodiment, waveform inversion is performed in the Laplace domain, but the waveform inversion is performed in the hybrid domain of the Laplace-Fourier domain using a complex damping constant.

For example, when the signal processor 103 processes data using waveform inversion in order to visualize subsurface velocity structure of a target region, the signal processor 103 first performs waveform inversion in the Laplace-Fourier domain to find a medium-wavelength or long-wavelength component for a velocity model of a target region, sets the velocity model as an initial velocity model for waveform inversion in the frequency domain, and performs frequency-domain waveform inversion to find a short-wavelength component for the velocity model.

FIG. 7 is a block diagram of the signal processor 103 according to an exemplary embodiment of the present invention.

A first input unit 201 is connected to the receiver 101 to receive time-domain data of the target region. The first input unit 201 transmits the received time-domain data to a transforming unit 202. For example, data which is transmitted from the first input unit 201 to the transforming unit 202 can be expressed by the above Equation 2.

The transforming unit 202 transforms the received time-domain data into Laplace-Fourier domain data. For example, the transforming unit 202 can transform the time-domain data into Laplace-Fourier-domain data using the complex damping constant as in Equation 1.

A modeled data generating unit 203 generates modeled data of the target region. For example, if parameters representing physical characteristics of the target region is input to a second input unit 204, a parameter storing unit 205 stores the parameters, and the modeled data generating unit 203 receives the parameters stored in the parameter storing unit 205, sets up an equation including the parameters, and solves the equation in the Laplace-Fourier domain to generate modeled data.

That is, an output of the transforming unit 202 means measured data, and an output of the modeled data generating unit 203 means estimated data.

A parameter update unit 206 compares an output (that is, measured data in the Laplace-Fourier domain of the transforming unit 202) of the transforming unit 202 with an output (that is, modeled data of the Laplace-Fourier domain) of the modeled data generating unit 203, thus updating the parameter stored in the parameter storing unit 205. At this time, the parameter update unit 206 can update the previously stored parameter to minimize a difference between the measured data and the modeled data.

The process will be described below in detail using equations.

The modeled data generating unit 203 can express a wave equation including parameters representing the characteristics of the target region as Matrix Equation 4, using a finite element method (FEM) or a finite difference method (FDM):

[Corrected under Rule 26 07.08.2009]
,(4)

where M denotes a mass matrix, C denotes a damping matrix, K denotes a stiffness matrix, u denotes a time-domain wavefield, f denotes a source vector, and · denotes differentiation with respect to time. Here, M, C and K are functions of the parameters representing physical characteristics of the target region.

If a solution u of Equation 4 is obtained, u can be modeled data of the target region. In order to obtain the solution in the Laplace-Fourier domain, Equation 4 can be transformed as follows:

[Corrected under Rule 26 07.08.2009]
,(5)

The solution of Equation 5 in the Laplace-Fourier domain can be simply obtained by factoring a complex impedance matrix S and performing a forward or backward substitution method or an iterative matrix solving method with respect to a source vector.

Then, the parameter update unit 206 compares u(s) obtained by the modeled data generating unit 203 with d(s) obtained by the transforming unit 202 to update the parameter stored in the parameter storing unit 205 so that the difference between u(s) and d(s) can be minimized. That is, the parameter can be updated in such a manner as to define an objective function E representing the difference between d(s) and u(s) and minimize the objective function E.

For example, a difference between d(s) and u(s) with respect to an i-th Laplace-Fourier frequency of a k-th receiver by a j-th source can be expressed using logarithms by Equation 6:

[Corrected under Rule 26 07.08.2009]
(6)

[Corrected under Rule 26 07.08.2009]
The objective function E can be defined by Equation 7: ,(7)

where nf denotes the number of frequencies for inversion, ns and nr denote the number of sources and receivers, respectively, and * denotes a complex conjugate.

Updating the parameter by minimizing the objective function E can be performed by a steepest descent method, and the steepest descent direction of an m-th model parameter pm for a given complex frequency can be expressed by Equation 8:

[Corrected under Rule 26 07.08.2009]
(8)

Equation 8 is normalized using a diagonal element of the pseudo-Hessian matrix which uses a virtual source vector, and parameters can be updated using the obtained steepest descent direction as follows:

[Corrected under Rule 26 07.08.2009]
(9)

where NRM denotes normalization, pml denotes an m-th model parameter which is once iterated, denotes a damping factor, and denotes a scale coefficient of a step length.

Model parameters which minimize the objective function E can be obtained by iteratively summing a scaled update increment of the steepest descent direction, starting from an initial parameter, as described above.

Since the parameters which minimize the objective function E can be regarded as identical to the physical characteristics of the target region through updating, an image data generating unit 207 can generate image data of the target region using the finally updated parameters.

Alternatively, the objective function E can be defined using a difference between the p-th power of d(s) and the p-th power of u(s), or using a difference of integrals of the p-th power of d(s) and the p-th power of u(s) with respect to p, wherein p can be smaller than 1.0. Also, any other numerical analysis techniques which minimize the objective function E can be used.

FIG. 8 is a flowchart of a subsurface velocity structure visualization method according to an exemplary embodiment of the present invention.

The subsurface velocity structure visualization method can be executed by the subsurface velocity structure imaging apparatus, or can be implemented in the form of a program stored on a recording medium which can be executed in a subsurface velocity structure visualization apparatus.

Referring to FIG. 8, in operation S101, measured data is received from the target region in step S101. The measured data may be seismic data propagated in the target region and can be expressed using a time-domain wavefield like Equation 2.

Then, in operation S202, the measured data is transformed to data in the Laplace-Fourier domain. The transformation process can be performed using Equation 1, and a complex number which includes a Laplace damping constant and angular frequency can be used as a damping constant s upon the transformation.

Then, in operation S103, a target region is modeled in step S103. For example, a matrix equation such as Equation 5 is set using parameters which represent characteristics of the target region, and Equation 5 is solved to generate modeled data. In operation S103, a process for receiving initial parameters for the target region can be additionally performed.

Subsequently, in operation S104, the objective function representing a difference between the measured data and the modeled data in the Laplace-Fourier domain is generated in step S104. For example, the objective function can be defined using a logarithmic value of each data like Equations 6 and 7. Alternatively, the objective function can be defined using the p-th power of each data or an integral of the p-th power of each data with respect to p.

Next, the initial parameters are iteratively updated to minimize the objective function in operations S105 and 106. Updating the parameters can be iterated until an objective function generated using the updated parameters is compared with and determined as equal to or less than a predetermined threshold value.

For example, if the objective function is determined as equal to or less than the predetermined threshold value in operation S105, operation S107 is performed, and if not, the parameters for generating the modeled data are updated in operation S106. At this time, the parameters can be updated to minimize the objective function.

When the parameters are finally updated through the above-described process, image data of the target region is generated using the parameters in operation S107. For example, when a velocity model of the target region is used as a parameter, the image data can be generated based on a velocity model which is finally obtained by iteratively updating the velocity model.

FIG. 9 shows a P-wave velocity model of a BP model according to an exemplary embodiment of the present invention.

Here, the BP model, which is disclosed in "The 2004 BP Velocity Benchmark" by F.J. Billette and S. Brandsberg Dahl, is a virtual subsurface velocity structure model which is used to inspect subsurface exploration technologies. That is, when it is assumed that the target region has a structure of FIG. 9, the effectiveness of the subsurface exploration technology can be tested by determining whether a structure which is identical or similar to that of FIG. 9 is obtained from measured data (this can be obtained from a data set of the BP model).

FIG. 10 shows a time-domain wavefield of the BP model.

FIG. 10 can be derived from an artificial data set of the BP model. For example, FIG. 10 corresponds to a seismogram obtained through field exploration of the target region.

FIG. 11 shows a wavefield in the Laplace-Fourier domain with respect to FIG. 10.

(a), (b) and (c) of FIG. 11 correspond to the cases where the Laplace damping constants are 0, 1 and 10, respectively. As the Laplace damping constant increases, the amplitude and phase spectrums of the Laplace-Fourier wavefield gradually smooth out.

FIG. 12 shows an example of an initial velocity model for a waveform inversion.

Referring to FIG. 12, the initial velocity model is set to have a linear velocity distribution according to depth.

(a), (b) and (c) of FIG. 13 show the steepest descent directions in the frequency domain, the Laplace domain, and the Laplace-Fourier domain, respectively.

In detail, (a) of FIG. 13 shows the steepest descent direction in the frequency domain, corresponding to short-wavelength information, (b) of FIG. 13 shows the steepest descent direction in the Laplace domain, corresponding to long-wavelength information, and (c) of FIG. 13 shows the steepest descent direction in the Laplace-Fourier domain, corresponding to medium-wavelength information which is a wavelength shorter than that of (b) of FIG. 13.

(a) of FIG. 14 shows a velocity model finally updated through waveform inversion in the frequency domain, (b) of FIG. 14 shows a velocity model finally updated through waveform inversion in the Laplace domain, and (c) of FIG. 14 shows a velocity model finally updated through waveform inversion in the Laplace-Fourier domain. The updated velocity model by the Laplace-Fourier-domain waveform inversion is the best among three results.

FIG. 15 shows another example of a velocity model finally updated through waveform inversion. The velocity model of FIG. 15 is a result obtained by performing waveform inversion once more in the frequency domain when the velocity model obtained in (c) of FIG. 14 is set as an initial velocity model.

As seen in FIGS. 14 and 15, the velocity model obtained according to the exemplary embodiments is very similar to that of FIG. 9 which is assumed as the true velocity model. Particularly, a subsurface velocity model similar to the true velocity model can be obtained by applying waveform inversion in the Laplace-Fourier