KR101290332B1 - seismic imaging apparatus utilizing macro-velocity model and method for the same - Google Patents

Abstract

Waves output from a transmitter source are related to seismic imaging techniques that image underground structures by processing received measurement data after propagating underground structures. The disclosed underground imaging method calculates the imaging data of the underground structure by waveform inversion using a macro-velocity model as an initial value. The long-wave velocity model, which is the initial value of the waveform inversion, calculates a velocity difference function, which is the difference between the actual velocity and the initial velocity model, as a function of the ratio of measured elastic wave scattering energy and modeling scattering energy based on the initial velocity model. And updating the initial velocity model to calculate the long wavelength velocity model.

Description

Seismic imaging apparatus utilizing macro-velocity model and method for the same

The present invention relates to a seismic imaging technique in which waves output from a source propagate underground structures and then process the received measurement data to image the underground structures.

A full waveform inversion technique is a method of studying underground structures using information recorded at the surface. However, it is very difficult to accurately recover a velocity model from measured data or even artificial data due to the presence of local minimums of the objective function or the lack of low frequency components in the measured data. Because of this, waveform inversion is a two-step process. In the first stage, the long wavelength velocity model is restored, and in the second stage, waveform inversion is performed to improve the resolution of the velocity model.

Several methods are known for obtaining long wavelength velocity models. The method by traveltime tomography disclosed in "Luo, Y. and Schuster, GT 1991. Wave-equation traveltime inversion.Geophysics 56, 645-653", or "Shin, C. and Cha, YH 2008. Waveform inversion in the Laplace domain. Geophysical Journal International 173, 922-931. 2008 ".

An object of the present invention is to reduce the amount of calculation in waveform inversion and to provide a clearer image.

Furthermore, an object of the present invention is to improve the resolution of a deep region image by improving the characteristics of the low frequency band.

Furthermore, an object of the present invention is to improve the convergence speed of the calculation by improving the initial velocity model in waveform inversion.

According to an aspect, a new method for restoring a macro or ultra-macro velocity model of an underground structure is proposed. This new aspect is based on the observation that reflected waves from fast structures have a large amplitude. Accordingly, it is assumed that the high speed structure of the underground structure can be restored using the energy calculated from the measured data. Furthermore, it is assumed that the energy of the seismogram observed is similar to the gravity anomaly in gravity inversion. Gravity shifts from each subterranean structure are similar to the energy of perturbed signals from each scatterer.

According to another aspect, wave field data is converted into energy data to recover high-velocity anomalies. Based on the Born theory described in "Lkelle, LT and Amundsen, L. 2005. Introduction to Petroleum Seismology. Society of exploration Geophysics. ISBN 1560801298", wave field data is scattered by underground point scatterers. It is assumed to be equal to the sum of the signals. A system of equations is built that relates the energy of the observed data to the energy of the scattered signals.

According to these and further aspects described below, the long-wave velocity model can be calculated only by calculation by analytical solution without the need for iterative calculation. In addition, by using the improved initial velocity model in waveform inversion, it is possible to avoid falling into the local minimum value, and the resolution of the deep underground image can be improved by improving the low frequency characteristic using the long wavelength velocity model.

1 is a flowchart illustrating a schematic configuration of an imaging method of an underground structure according to an exemplary embodiment.
2 is a block diagram illustrating a schematic configuration of an imaging apparatus of an underground structure according to an exemplary embodiment.
3A and 3B are graphs showing the values of three different velocity disturbance functions Δv of different point-scattering circles.
4A and 4B are graphs showing the peaks of the scattered wave amplitude versus the magnitude of the disturbance rate of the signals measured at different source and receiver pairs.
5 illustratively shows the propagation path of a single point-scattered signal.

The foregoing and further aspects of the present invention will become more apparent through the following embodiments. In the following, these aspects will be described in detail with reference to the embodiments described with reference to the accompanying drawings.

1 is a flowchart illustrating a schematic configuration of an imaging method of an underground structure according to an exemplary embodiment. In the underground imaging method according to an embodiment, imaging data of the underground structure is calculated by waveform inversion using a macro-velocity model as an initial value. According to one aspect, a step of calculating a speed difference function that is a difference between an actual speed and an initial speed model from a function of a ratio of measured elastic wave scattering energy and modeling scattering energy based on the initial velocity model is performed. S10 to S40, and calculating a long wavelength velocity model by updating an initial velocity model using the velocity difference function (S60).

First, “Mitra, S.K. 2005. Digital signal processing: A computer-based approach. McGraw-Hill College. According to the definition proposed by Mitra in ISBN 0073048372 ”, the energy of the signal x (t) is defined as zero-lag autocorrelation.

(1)

(One)

According to the Lippmann-Schwinger integration described in Lkelle and Amundsen, 2005, cited above, the actual medium can be divided into two parts. One is a background velocity model and the other is a perturbed model, which is a scattered model. Waves propagating directly from a source to a receiver are called direct waves or unperturbed waves, and waves propagating through scattering sources are scattered waves or perturbed waves. It is called a wave. For models with point scatterers, wave field data can be separated by the Lippmann-Schwinger equation as follows:

(2)

(2)

Where u (t) is the acquired wave field using one source-receiver pair, u _i (t) is the non-random wave component of wave field u (t), and u _p (t) is the wave containing scattering wave component It is an intestinal ingredient.

Scattering can occur at any point underground and u _p (t) can be approximated by the sum of scattering waves by a single point scattering source, according to the Born approximation equation. Thus u _p (t) can be written as

(3)

(3)

Where s (x, t) is the 1st-order point-scattered signal and is the position of the point scatterer in the x domain Ω.

According to one aspect, u _i (t) includes signals preceding the bottom arrival, as well as direct arrivals, as well as seafloor reflections, such as scattered signals from below the bottom. Are assumed not to be included. Therefore, you can remove them from the data set.

Where s ( x , t) is an approximation because the computational load is severe. To define an appropriate approximation function, the sensitivity of the amplitude of the wave field to velocity disturbances using a point scattering source is empirically observed. 3A and 3B are signals obtained by using a fourth order finite difference algorithm as a function Δv of velocity disturbance. As shown in FIG. 3B, the magnitude of the velocity disturbance does not affect the shape of the scattered signal but only the magnitude. 3A shows the peak amplitude value of the scattering wave as a function of velocity disturbance Δv for five scattering sources with different source and receiver positions. Normalizing each data to remove the effect of geometrical spreading is shown in FIG. 3B. The resulting peak amplitude R [Δv] is a function of only Δv. From these two simple analyzes, we conclude that the amplitude variation of a single point scatterer signal can be divided into geometric variance and velocity disturbance components. The effect of geometric dispersion is “Kuvshinov, BN and Mulder, WA 2006. The exact solution of the time-harmonic wave equation for a linear velocity profile. Geophysical Journal International 167, 659-662. ”. According to this observation, S (x, t) can be expressed as follows.

(4)

(4)

Where x, x _s and x _r are the positions of the point scatterer, the transmitter and the receiver, respectively.

λ (x _b ; x _a ) is the ratio of amplitude attenuation due to geometric dispersion during wave progression from x _a to x _b , and Δv (x) is the difference between the actual velocity and the background velocity model at the scattering source position x. R [Δv (x)] is a function of relating the rate of energy scattering to the rate difference. λ (x _b ; x _a ) is calculated by the analysis solution described later. δ (tt _x ) is a delta function. T (x, x _s ) and T (x, x _r ) are the traveltimes of the scattered signal from point x to source x _s and point x to source x _r , respectively. A ₀ is the initial amplitude value of the transmission source.

However, delta function signal approximation is not suitable for expressing the energy of the scattering signal. Although the amplitude of the delta function is equal to the maximum amplitude of the disturbed signal, it can be easily seen that the energy of the scattering signal of equation (1) cannot be described by the delta function assumption. Therefore, the equation (4) is modified as follows by changing the approximation of the first-order scattering signal from the delta function to the box-shape signal.

(5)

(5)

This signal exists only between t _x and t _x + ε. By properly choosing the value of ε, one can assume that s (x, t) has the same energy as s (x, t) of a single point-scattering signal.

Equation (3) is redefined using pseudo primary point scattering signal s (x, t).

(6)

(6)

5 illustratively shows the propagation path of a single point-scattered signal. During propagation, the amplitude of the signal decreases in proportion to the propagation distance. Λ was calculated using a linearly increasing velocity model as background velocity and using the amplitude of the analytical solution according to the method disclosed in Kuvshinov and Mulder, 2006. In addition, this solution is not a halfspace solution but is used as an approximation. The interpretive solution can be expressed as

(7)

(7)

here,

(8)

(8)

to be.

Where v ₀ is the velocity at the surface, z _max is the maximum depth, v _max is the bottom velocity of the linear increasing velocity model, and x _a (x _s , z _s ) and x _b (x, z) are the source and receiver Are the position vectors of and r is the distance between the transmitter and the receiver period.

Equation (6) shows that u _p (t) can be approximated by the sum of the point-scattered wave fields. From the definition of the pseudo-first point-scattered signal and the energy of equation (1), the equilibrium equation of the scattered energy in each transmission cost receiver can be derived as follows.

(9)

(9)

Where e ^{s , r} is the energy of the signal recorded at position r for source s. However, it is difficult to calculate the system matrix equation of equation (9), and accurate traveltime is required. To solve this problem, it is assumed that the total energy can be approximated by the sum of the energy of each scattered signal. For this assumption to be satisfied, the interference between each scattered signal must not affect the energy distribution of the data. Sparse parameterization and high frequency approximation were used to minimize interference between the scattered signals by plotting them as box-shaped signals with short durations. Accordingly, equation (9) can be summarized by the following equation (10).

(10)

(10)

According to one aspect, the calculation of the velocity difference function may include approximating the measured or observed seismic signal to the sum of the scattered waves to calculate the measured seismic scattering energy (step S10 of FIG. 1). That is, on the right side of equation (10), energy is calculated at each receiver location for one shot and can be calculated as follows from the observed wave field data.

(11)

(11)

Here it was denoted by modifying the u _p (t) for each source and sujingi to ^{_{u p s, r (t)}} . Now we can calculate n _r × n _s scattered energy data. The left side of equation (10) is derived as follows. The energy from the scattered signal is expressed as follows using the definition of energy.

(12)

(12)

Where time _x is located at the center of the spherical signal. In addition, it is assumed that the amplitude varies depending on R [? ( X )] and the propagation distance in Equation (4). In this case, the sum of scattering energy is expressed as follows.

(13)

(13)

Finally, the equilibrium equation of scattering energy is summarized as follows.

(14)

(14)

Dividing Equation (14) by 2 ◎ A ₀ ² gives us an interesting characteristic:

(15)

(15)

Here, the proper choice of ε is important for the square signal assumption to be justified. The value of ε depends on the source signature. If the source characteristic is invariant for all blasts, then the value of epsilon is constant. Therefore, ε can be ignored in this algorithm which determines the relative proportion of the speed difference.

According to one aspect, the calculation of the velocity difference function may comprise the step of calculating the energy of the pseudo-scattered wave signals from a simple velocity model, for example a linear increasing velocity model (step S20 of FIG. 1). In other words, the integral value of the square of the product of λ (x; x _S ) and λ (x; x _r ) at the left side in equation (15) can be calculated by an analysis solution of equation (7).

According to another aspect, the calculation of the velocity difference function may include calculating a relative velocity shift ratio that is a ratio of the measured value of the seismic scattering energy and the energy value of the pseudo-scattering wave signals (step S30 of FIG. 1). have. That is, equation (15) can be expressed as follows.

(16)

(16)

Where K is an (n _s × n _r ) × (n _x × n _z ) matrix, r is a vector with (n _x × n _z ) components, and e is a vector with (n _s × n _r ) components N _s and n _r are the number of transmitters and receivers, respectively. This system equation can be solved using regular equations as

(17)

(17)

r can be solved using only diagonal components, taking into account computational resource limitations. By using only the diagonal components, the algorithm is very fast and can be processed at low cost. We can then see where the speed shift is located by calculating the relative ratio, r . However, since we used only the diagonal components in the matrix, the negative value of r cannot be obtained.

According to one aspect, the calculating process of the speed difference function includes calculating the speed difference function using the relationship between the relative speed shift ratio and the speed difference function shown in step S40 of FIG. 1. The velocity shift ratio Δv ( x ) can be solved with reference to the equation r (x) = R [Δv ( x )] ² and the graph of FIG. 4B. The computed shift r ( x ) shows where the high speed structure is located. However, the results may be unclear and widespread. The maximum position in the resulting value may be exactly the position of the structure, but the surrounding area including the fast structure may also have a large value. Based on this observation, it seems desirable to abbreviate the range of results through scaling.

According to an aspect, the seismic imaging method may further include normalizing the speed difference function after calculating the speed difference function (step S50 of FIG. 1). The size of the anomalous body can be controlled by the following equation.

(18)

(18)

Where N | · | Is a normalization operation of a value, and q is a scaling constant. As the q value increases, the size of the shift volume decreases. In this example, q = 2 was used. Through this process, more reasonable result of speed update can be obtained.

According to one aspect, the calculating process of the speed difference function includes calculating the speed difference function using the relationship between the relative speed shift and the speed difference function (step S60). In this embodiment, the velocity model is calculated by updating the initial velocity, ie, the simple linear increasing velocity model, with the scaled shift value as follows.

(19)

(19)

Where V _r (x) is the speed calculated according to the invention, V ₀ (x) is the initial speed and γ is the step length of the update. In order to properly determine the values of q and γ, several linear searches were performed to calculate the difference between seismograms obtained using actual and calculated velocity models.

2 is a block diagram illustrating a schematic configuration of an imaging apparatus of an underground structure according to an exemplary embodiment. As shown, an underground imaging apparatus according to an embodiment calculates a velocity difference function that is a difference between an actual velocity and an initial velocity model from a function of a ratio of measured acoustic wave scattering energy and modeling scattering energy based on the initial velocity model. Based on the velocity difference function calculation unit 100, the velocity model calculation unit 500 for updating the initial velocity model according to the velocity difference function, and calculating a long wavelength velocity model, and an underground structure by waveform inversion based on the calculated long wavelength velocity model. It may include an underground structure imaging unit 700 for calculating the imaging data of the.

The underground structure imaging unit 700 continuously updates the modeling parameters in a direction of minimizing the residual function between the modeling data and the measured data, which are solutions of the wave equation, and calculates the modeling parameters which are speed or density values of the underground structure.

According to another aspect, the velocity difference function calculating unit 100 approximates the observed seismic signal to the sum of the scattering waves and calculates the observed seismic energy, which is the scattering energy, with the observed seismic energy calculating unit 130. The modeling energy calculator 110 calculates the energy of the pseudo-scattered signals from a simple velocity model, the observed seismic energy value and the pseudo-scattered wave energy value. Scattered wave energy ratio calculation unit 150 for calculating a relative velocity anomalies ratio that is a ratio of, and an output unit 170 for calculating a velocity difference function using the relationship between the relative velocity deviation ratio and the speed difference function. It may include.

As described above, the observed acoustic wave energy calculator 130 calculates the observed acoustic wave energy based on Equation (11). In addition, the modeling energy calculator 110 calculates the energy of the pseudo-scattered wave from the simple velocity model, for example, the linear increase velocity model. The modeling energy calculation unit 110 calculates the energy of the pseudo-scattered wave represented by the K matrix in Equation (16) using the analysis solution of Equation (7). Subsequently, the scattered wave energy ratio calculator 150 calculates a relative speed shift ratio that is a ratio between the observed acoustic wave energy value and the pseudo-scattered wave energy value based on Equation (17). Thereafter, the output unit 170 calculates the speed difference function using the relationship between the relative speed shift ratio and the speed difference function illustrated in FIG. 4A.

According to another aspect, the imaging apparatus of the underground structure may further include a normalizer 300 that normalizes the speed difference function calculated by the speed difference function calculator 100 and outputs the normalized speed difference function to the speed model calculator. The normalization unit 300 normalizes the speed difference function based on, for example, Equation (18).

The present invention has been described above with reference to the embodiments described with reference to the accompanying drawings, but is not limited thereto and is construed to cover many modifications that can be obviously derived therefrom. The appended claims are intended to cover such modifications.

100: speed difference function calculation unit 110: modeling energy calculation unit
130: observed seismic energy calculation unit 150: scattering wave energy ratio calculation unit
170: output unit 300: normalization unit
500: speed model calculator 700: underground structure imaging unit

Claims (6)

Underground structure imaging method that uses the macro-velocity model as an initial value and calculates the imaging data of the underground structure by waveform inversion. The long-wave velocity model, which is the initial value of the waveform inversion, is used. This:
Calculating a velocity difference function that is a difference between the actual velocity and the initial velocity model from a function of the ratio of the measured acoustic wave scattering energy and the modeling scattering energy based on the initial velocity model;
Calculating a long wavelength velocity model by updating an initial velocity model with the velocity difference function;
Imaging method of the underground structure calculated through a process comprising a.
The method of claim 1, wherein calculating the speed difference function comprises:
Approximating the observed seismic signal to the sum of scattered waves to calculate the observed seismic energy, the scattering energy;
Calculating the energy of the pseudo-scattered signals from a simple velocity model;
Calculating a relative velocity anomalies ratio, which is a ratio of the value of the observed seismic energy to that of the quasi-scattered wave;
Calculating a speed difference function using the relationship between the relative speed shift ratio and the speed difference function;
Imaging method of the underground structure comprising a.
The method of claim 1, wherein after calculating the speed difference function:
Normalizing the calculated speed difference function;
Imaging method of the underground structure further comprising.
In the imaging apparatus (seismic imaging apparatus) of the underground structure,
A velocity difference function calculation unit that calculates a velocity difference function that is a difference between the actual velocity and the initial velocity model from a function of the ratio of the measured acoustic wave scattering energy and the modeling scattering energy based on the initial velocity model;
A speed model calculator for updating the initial speed model according to the speed difference function to calculate a long wavelength speed model;
An underground structure imaging unit configured to calculate imaging data of the underground structure by waveform inversion based on the calculated long wavelength velocity model;
Imaging device of the underground structure comprising a.
The method of claim 4, wherein the speed difference function calculating unit:
An observed seismic energy calculator for approximating the observed seismic signals to the sum of the scattering waves to calculate the observed seismic energy;
A modeling energy calculator for calculating the energy of the pseudo-scattered signals from a simple velocity model;
A scattering wave energy ratio calculating unit for calculating a relative velocity anomalies ratio which is a ratio of the observed acoustic wave energy value and the energy value of the pseudo-scattering wave;
An output unit for calculating a speed difference function using a relationship between a relative speed shift ratio and a speed difference function;
Imaging device of the underground structure comprising a.
5. The method of claim 4,
A normalizer for normalizing the speed difference function calculated by the speed difference function calculator and outputting the normalized speed difference function to the speed model calculator;
Imaging apparatus of the underground structure further comprising.

Images