WO2016060513A1 - 탄성파 다성분 자료에 대한 중합전 egs 구조보정 방법 - Google Patents

탄성파 다성분 자료에 대한 중합전 egs 구조보정 방법 Download PDF

Info

Publication number
WO2016060513A1
WO2016060513A1 PCT/KR2015/010951 KR2015010951W WO2016060513A1 WO 2016060513 A1 WO2016060513 A1 WO 2016060513A1 KR 2015010951 W KR2015010951 W KR 2015010951W WO 2016060513 A1 WO2016060513 A1 WO 2016060513A1
Authority
WO
WIPO (PCT)
Prior art keywords
wave
egs
propagator
frequency
correction
Prior art date
Application number
PCT/KR2015/010951
Other languages
English (en)
French (fr)
Inventor
김병엽
변중무
설순지
이호영
Original Assignee
한국지질자원연구원
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 한국지질자원연구원 filed Critical 한국지질자원연구원
Priority to US15/508,957 priority Critical patent/US20170299745A1/en
Publication of WO2016060513A1 publication Critical patent/WO2016060513A1/ko

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/362—Effecting static or dynamic corrections; Stacking
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S17/00—Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/282—Application of seismic models, synthetic seismograms
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/284—Application of the shear wave component and/or several components of the seismic signal
    • G01V1/286—Mode conversion
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/32—Transforming one recording into another or one representation into another
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V9/00—Prospecting or detecting by methods not provided for in groups G01V1/00 - G01V8/00
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/05—Geographic models
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10K—SOUND-PRODUCING DEVICES; METHODS OR DEVICES FOR PROTECTING AGAINST, OR FOR DAMPING, NOISE OR OTHER ACOUSTIC WAVES IN GENERAL; ACOUSTICS NOT OTHERWISE PROVIDED FOR
    • G10K15/00—Acoustics not otherwise provided for
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/50—Corrections or adjustments related to wave propagation
    • G01V2210/51—Migration
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/50—Corrections or adjustments related to wave propagation
    • G01V2210/51—Migration
    • G01V2210/512—Pre-stack

Definitions

  • the present invention relates to a method for migrating an acoustic wave, and more particularly, to a conventional scalar generalized-screen (SGS) capable of quickly calculating wave propagation in a medium having a horizontal velocity change.
  • SGS generalized-screen
  • an elastic generalized membrane can be used to efficiently represent the behavior of an elastic wave that propagates the interface of the underground medium and undergoes a mode change between P and S waves.
  • Screen (EGS) One-way wave equation using wave propagator relates to a method for prestack depth migration before polymerization.
  • the present invention is to extend the Taylor series expansion of the vertical slowness term of the wave propagator to the second order to calculate a higher accuracy wave field in a medium of complex structure A method for correcting pre-polymerized EGS structure for data.
  • the present invention unlike the conventional techniques used after first separating the input data into the P wave and the S wave wave field when the structure correction of the multi-component data structure, by including the mode separation operator in the propagator,
  • the present invention relates to a pre-polymerization EGS structure correction method for an elastic wave multi-component data configured to generate P wave and S wave image cross sections by using a shot gather as a structure correction input without having to separate the S wave.
  • the present invention because the polarity is reversed at the reflection point in the case of the seismic S-wave sound source collection (shot gather), if used as it is in the structure correction, the continuity of the event is lowered and the quality of the final cross section is deteriorated
  • a module for correcting the polarity switching in the frequency-frequency domain before S-wave imaging it is possible to improve the quality of the S-wave structure correction image.
  • seismic structure migration involves moving all the primary reflection events on the stack section into place, and by collapsing the diffraction curves to increase spatial resolution and obtain images of underground structures. It means the process to make it.
  • the correction position is determined by the travel time and the velocity of the reflected wave, and the structure correction method is divided largely into the Kirchhoff structure correction based on the integral solution of the wave equation, and the finite difference ( It is divided into structural correction using finite difference method and frequency-wave frequency (FK) structural correction using frequency and wave number.
  • the structure correction method is divided largely into the Kirchhoff structure correction based on the integral solution of the wave equation, and the finite difference ( It is divided into structural correction using finite difference method and frequency-wave frequency (FK) structural correction using frequency and wave number.
  • Kirchhoff's structural correction method calculates an appropriate diffraction hyperbola at each point on the polymerized cross section based on the dashed line theory and verifies the diffraction hyperbola at that point.
  • the sum of all the sample values in the sample is taken as the amplitude of the corresponding point on the structural correction cross section.
  • the calculation speed is fast and relatively accurate results are used in the industry.However, for complex structures, the dashed line theory is not applied properly. There is a drawback to failure.
  • the finite-difference structural correction technique is generally the most widely used reverse time migration method using the full two-way wave equation, which is relatively accurate even in complex structures. It has the advantage of deriving, but it has the disadvantage of taking too much calculation time.
  • the structural correction technique using the one-way wave equation has the disadvantage of the Kirchhoff method, which may fail to image in a complex structure as described above, and the inverse time that requires a lot of calculation time.
  • the proposed method is designed to consider only the upgoing wave propagating in one direction by approximating the bidirectional wave equation, so that it can be solved by finite difference in the time domain or the frequency-frequency (FK) domain. It is a way of spreading waves.
  • the one-way wave equation is used as a propagator for propagating waves, and the final structural correction section is obtained by applying the wave field generated in the sound source region and the wave field generated in the receiver region to the imaging conditions. Equation restructuring technique.
  • the wave equation can be processed, but the earth's underground medium is strictly water-acoustic. It is not a medium but an elastic media with complex and heterogeneous properties.
  • the acoustic wave returned from the elastic medium must be detected with a three-component geophone that can detect the horizontal and vertical displacement (or velocity, acceleration) of the medium to accurately describe its behavior.
  • the data should be processed using the elastic wave equation.
  • the structure correction technique using the elastic wave equation is similar to the structure correction technique using the scalar wave equation, but multicomponent seismic data is used as the input wave field, and the wave equation is also the elastic rather than the acoustic wave equation.
  • the mode transition the transition of P-to-S-wave or S-to-P-wave mode
  • the attenuation of waves ie, P-wave and S-wave
  • the conventional EGS propagator whose vertical slow term is calculated only up to the first term is extended to the second term, so that even more accurate waves can be obtained even in a complex model or a model with a high speed change horizontally. It is desirable to provide a new structural correction algorithm that can be configured to correct the polarization reversal of S-waves by implementing the propagation of the E-wave. I have not been given a way.
  • the present invention seeks to solve the problems of the prior art as described above, and therefore an object of the present invention is to solve the problem of the prior art EGS propagator, in which the vertical slowness operator is limited to approximate only the first term.
  • the vertical slow term In order to improve the performance of EGS radio wave propagators, we extend the vertical slow term to the second term to realize more accurate wave propagation even in complex models or horizontally varying speed models. It is to provide a method for correcting EGS structure before polymerization.
  • EGS propagator Including internally implemented PS separation module, P- and S-wave image sections can be generated by using shot gather as a structure correction input without the need to separate the input multi-component wave field into P and S waves. It is to provide a method of pre-polymerization EGS structure correction for the composed seismic multicomponent data.
  • another object of the present invention is that, in the case of the seismic S-wave sound source collection (shot gather), since the polarity is reversed at the reflection point, if it is used for structural correction as it is, the continuity of the event is lowered.
  • a module for correcting polarity switching in the frequency-wavelength region before S wave imaging is added to improve the quality of the S structure correction image.
  • the purpose of this paper is to provide a pre-polymerization EGS structural correction method for multi-component seismic data.
  • an EGS (Elastic Generalized-Screen) wave propagator for expressing the behavior of a seismic wave through the mode switching between P and S waves while propagating the interface of the underground medium
  • the pre-polymerization EGS structure correction method for the acoustic wave multi-component data configured to generate P- and S-wave image cross sections, the acoustic wave multi-component data to be analyzed are input to the source and receiver.
  • a preprocessing step of establishing a model and determining a frequency band to be calculated through a Fourier transform A sound source region processing step of calculating a wave propagation from source in a sound source region for each frequency band by using the EGS radio wave; A receiver region processing step of calculating a backward propagation from receiver in a receiver region for each frequency band using the EGS propagator; A structure correction step of integrating the results calculated in the sound source region processing step and the receiver region processing step through cross correlation and performing a structural correction by an imaging condition; And an output step of outputting the structure-corrected image data in the structure correction step.
  • the EGS radio wave is characterized in that it is configured to use the EGS radio wave shown in the following equation.
  • the sound source region processing step may further include: separating a sound source wavefield by a mode separation operator ([M 0 ] ⁇ 1 ) in the EGS propagator; Calculating screen and mode coupling in the frequency-space (fx) region; Performing a Fourier transform on the spatial variable (x) in the calculated result; Calculating an extrapolated wave field in the frequency-frequency region fk; Storing each mode for migration, and recomposing the sound source wave field by a mode combining operator (M 0 ) in the EGS propagator; And performing an inverse Fourier transform on the reconstructed sound source wave field.
  • a mode separation operator [M 0 ] ⁇ 1 ) in the EGS propagator
  • the receiver region processing step separating the receiver wavefield by a mode separation operator ([M 0 ] -1 ) in the EGS propagator; Calculating screen and mode coupling in the frequency-space (fx) region; Performing a Fourier transform on the spatial variable (x) in the calculated result; Calculating an extrapolated wave field in the frequency-frequency region fk; Storing each mode for migration and recomposing the oscillator wave field by a mode combining operator (M 0 ) in the EGS propagator; And performing an inverse Fourier transform on the reconstructed waveguide wave field.
  • a mode separation operator [M 0 ] -1
  • M 0 mode combining operator
  • the sound source region processing step and the receiver region processing step may be configured to allow parallel processing by allocating processing for each frequency band using a plurality of processors, thereby reducing the overall processing time. It is characterized in that it is configured to.
  • the structure correction step is characterized in that it is configured to use the imaging condition shown in the following equation as the imaging condition.
  • I ij is the final imaged image
  • Is a scalar source wavefield in the Fourier domain
  • ⁇ I and j are the vector wavefields of the sound source and the oscillator, and mean the horizontal component (x-component) and vertical component (z-component) in the multi-component seismic data, respectively.
  • the method includes a polarity correction step of correcting the S-wave polarity reversal phenomenon in which the polarity is changed at the reflection point of the medium boundary by obtaining the reflection angle at the reflection point in the frequency-wavelength region using the following equation. It is characterized in that it further comprises.
  • the shot gather is directly input to the structure correction input without the need to separate the multicomponent data into P and S waves. It is configured to generate P- and S-wave image cross sections, and to improve the S-wave structure correction image quality by performing correction for polarity switching in the frequency-frequency domain before S-wave imaging.
  • a pre-polymerization EGS structural correction system is provided for characterized seismic multicomponent data.
  • the vertical slowness operator is extended to the second term to realize more accurate wave propagation even in a complex model or a horizontally variable speed model, thereby further improving the performance of the EGS propagator.
  • a P wave and a shot wave is directly used as a structure correction input without having to separate the input multicomponent wave field into P and S waves, including a PS separation module implemented inside the EGS propagator.
  • EGS structure correction method is provided for pre-polymerization of ESA multi-component data configured to generate S-wave image cross-sections. Computation and structure are complicated by separating wave field into P wave and S wave as input before structure correction. It is possible to solve the problems of the conventional structure correction techniques, which had disadvantages.
  • the present invention before the S-wave imaging, by adding a module for correcting the polarity switching in the frequency-wavelength region to improve the quality of the S-wave structure correction images
  • polarity is reversed at the reflection point in the case of seismic S-wave collection, so if it is used for structure correction as it is, the continuity of the event is reduced and the quality of the final cross section is reduced. It is possible to solve the problems of the structural correction techniques of the prior art that had a problem of deterioration.
  • 1 is a diagram illustrating the concept of an EGS radio wave.
  • FIG. 2 is a diagram showing a wave field propagated to an EGS propagator according to an embodiment of the present invention in a zero-pertubation medium.
  • FIG 3 is a view showing a wave field propagated to an EGS propagator according to an embodiment of the present invention in a medium in which microvariation exists.
  • FIG. 5 is a diagram illustrating the concept of an EGS structure correction algorithm using the EGS propagator shown in FIG. 1.
  • FIG. 6 is a flowchart schematically illustrating an overall configuration of an EGS structure correction algorithm according to an embodiment of the present invention implemented based on the conceptual diagram illustrated in FIG. 5.
  • FIG. 7 is a flowchart schematically illustrating an overall configuration of a process of implementing MPI for parallel processing of an EGS structure correction algorithm according to an embodiment of the present invention shown in FIG. 6.
  • FIG. 8 is a diagram illustrating a comparison of structural correction results of applying and not applying a polarity inversion correction method using imaging conditions of an EGS structure correction algorithm according to an embodiment of the present invention.
  • FIG. 9 is a diagram illustrating a P wave velocity structure of a previously generated two-dimensional cross-sectional model for verification of an EGS structure correction algorithm according to an exemplary embodiment of the present invention.
  • FIG. 10 is a view showing images of final PP and PS cross sections derived by applying the EGS structure correction technique according to an embodiment of the present invention.
  • the present invention in order to solve the problem of the EGS propagator of the prior art, in which the vertical slowness operator is approximated only to the first term, as described below, the present invention extends the vertical slow term to the second term, which is complicated.
  • the present invention relates to a pre-polymerization EGS structure correction method for multi-stage seismic data, which is configured to realize more accurate wave propagation even in a model or a horizontally varying speed model.
  • the present invention as described below, to solve the problems of the structural correction techniques of the prior art had a disadvantage in that the calculation and structure is complicated by separating the wave field into the P wave and the S wave before use as an input, Including the PS separation module implemented inside the EGS propagator, P- and S-wave image sections can be generated by using shot gather as a structure correction input without the need to separate the input multi-component wave field into P and S waves.
  • a method for correcting pre-polymerized EGS structure for elastic wave multicomponent data was described below, to solve the problems of the structural correction techniques of the prior art had a disadvantage in that the calculation and structure is complicated by separating the wave field into the P wave and the S wave before use as an input, Including the PS separation module implemented inside the EGS propagator, P- and S-wave image sections can be generated by using shot gather as a structure correction input without the need to separate the input multi-component wave field into P and S waves.
  • the present invention in the case of the seismic S-wave sound source collection (shot gather), due to the phenomenon that the polarity is reversed at the reflection point, if used as it is for structural correction, the continuity of the event is reduced and the final cross section
  • a module for correcting the polarity shift in the frequency-wavelength region before the S-wave imaging is added to improve the quality of the S-wave structure correction image.
  • a method for correcting pre-polymerized EGS structure for multi-stage seismic data that is configured to improve.
  • a propagator is an operator or function such as an engine used to numerically set the earth model and then numerically simulate the progress of the seismic waves on the model. Therefore, the better the propagator is, the better it can be described in numerical modeling that the actual seismic wave propagates inside the earth.
  • the propagator presented in 2003 is a vertical slow term (the function used when the actual wave propagates downwards) calculated only up to the first order, which is, in other words, an incorrect propagator,
  • the propagator presented in the above-mentioned prior art document is merely mathematically deriving the propagator and has not been implemented by the structure correction algorithm.
  • the present inventors improve the conventional propagator calculated up to the first order to implement the structure correction algorithm, and at the same time, the EGS propagator is extended to the secondary order for more accurate wave propagation, thereby numerically improving the performance of the propagator.
  • a structure correction algorithm was implemented using the improved radio wave, and at this time, the P wave and the S wave were separated from the radio wave so that it could be used for structure correction and a module for correcting the polarity reversal of the S wave.
  • the performance of the EGS propagator can be further improved.
  • Le Rousseau 2001 (see Ref. 9), uses the Fourier domain split step (Stoffa et al., 1990, see Ref. 10) technique and phase-screen method (Wu and Huang, 1992, see reference 11).
  • the technique has been generalized to allow more accurate propagation in horizontally varying speed media, and this method has been termed a generalized-screen propagator, Le Rousseau and De Hoop, 2003 ( Reference 8 expands this into an elastic wave.
  • Equation 1 An elastic thin-slab propagator exhibiting an upward wave propagating upward and a downward wave propagating downward between thin slabs is represented by Equation 1 below. It is expressed in the form of pseudo-differential operator ( ⁇ DO).
  • Equation 1 expressed as a pseudo differential operator is an operator that extrapolates the wave field in the depth direction by shifting the phase using the vertical slow term ( ⁇ ) of the medium in the frequency-wave region.
  • x the Fourier variable ( ⁇ v ) are expressed as a single term.
  • the subscript ⁇ is a right symbol of a pseudo differential operator.
  • the separation of the spatial and frequency domain parameters of the elastic unidirectional thin film propagator by approximation is the core of the elastic generalized-screen propagator considering the horizontal velocity change.
  • Equation 3 For the characteristic operator A: Has the relationship as shown in
  • Equation 4 a symbol of A
  • the subscript indicates the symbol order and the superscript indicates the order of medium contrast.
  • a principal symbol matrix (a [2]) can be calculated.
  • ⁇ 22 , ⁇ 12 , ⁇ 13 and ⁇ 32 can be obtained by exchanging ⁇ 1 and ⁇ 2 from ⁇ 11 , ⁇ 21 , ⁇ 23 and ⁇ 31 , respectively, and use the Taylor series to the third term and k s -1, k s -1/2, ⁇ 1/2, -1/2 and expand ⁇ (expanding) to replace them by the aforementioned formula to rearrangement against the epsilon ( ⁇ ) the order (order) By (rearrange), a symbol of higher order can be calculated.
  • the microvariation of can be expressed as Equation 6 below.
  • Equations 7 below are terms independent of epsilon epsilon in Equation 5, and thus include the background term of the propagator.
  • ⁇ ⁇ , 22 , ⁇ ⁇ , 12 , ⁇ ⁇ , 13 and ⁇ ⁇ , 32 represent ⁇ 1 and ⁇ 2 , respectively, such as ⁇ ⁇ , 11 , ⁇ ⁇ , 21 , ⁇ ⁇ , 23 And ⁇ ⁇ , 31 .
  • Equation 8 is rearranged by the first epsilon oredr term after Taylor development of Equation 5, and Equation 9 is expressed by Equation 5 It was rearranged by the second epsilon term after Taylor development, where certain terms corrected errors in some of Le Rousseau's deviations (see Le Rousseau's, 2001, Ref. 9).
  • M 0 serves to couple the P and S waves separated after thin film propagation.
  • Equation 11 and Equation 12 are substituted into Equation 10, and Equation 8 By using, as shown in Equation 14 below 1st expansion of Can be obtained.
  • Diagonal components of the matrix are involved in the propagation of P and S waves, and the remaining components (off-diagonal entries) are involved in mode coupling, and all calculations are performed without separating the spatial variable ( ⁇ ) in the Fourier domain. There is a difficulty in performing the Fourier transform at the node to move between the space and the frequency domain.
  • the EGS propagator does not perform the Fourier transform by the number of nodes, but converts only the number of terms. It can be done quickly.
  • lambda represents the number of terms in [Equation 14] and [Equation 15].
  • Is calculated in the frequency-space (fx) region as the screen term Is computed in the frequency-frequency region fk with a phase-shift term.
  • N is a normalizing operator used to prevent inaccurate amplification or attenuation of energy due to wavenumber or propagation angle due to truncation errors during indefinite expansion (see Le Rousseau and De Hoop, 2003). See Document 8.
  • FIG. 1 is a diagram illustrating the concept of an EGS radio wave, and conceptually illustrates the calculation of Equation 17 described above.
  • the wave fields separated into P and S waves have two background media (the background velocities of the P and S waves, respectively).
  • Phosphorus Medium Thin film of phosphorous medium).
  • Second wave of S wave In the case of a medium, Not The use of is because the branch point may occur because the speed of the S wave is lower than the background speed of the P wave.
  • the P wave is The medium with the background velocity of ) Only contributes to mode switching and is not actually used in the next thin film.
  • the S wave Propagates a medium with a background velocity of. ) Is used only for mode switching and is not used for the next act.
  • P wave propagating medium with background speed of The S wave propagating through the medium with the background velocity becomes P wave and S wave actually proceeding to the next film.
  • the P wave and the S wave wave fields passing through the thin film are combined by M 0 to become the original wave fields Vx and Vz, and in the space region (FX region). To calculate the screen term.
  • the P wave and the S wave stored in each depth are stored as an array of wave depths in each array, and are used as input of imaging conditions in the structure correction.
  • this process becomes a key step in which multi-component seismic data can be used as input for structural correction, unlike other structural correction techniques in which the input wave field must be separated into P and S waves before structural correction.
  • FIG. 2 is a view showing a wave field propagated to an EGS propagator according to an embodiment of the present invention in a zero-pertubation medium.
  • the inventors of the present invention in order to verify the impulse response (impulse response) of the EGS propagator according to the embodiment of the present invention implemented by the algorithm shown in [Equation 17] and Figure 1, the P wave speed is 2,100 Impulse response after 0.2 s by generating a vertical point source at the surface on an isotropic homogeneous model with m / s, S wave velocity of 1,050 m / s, and density of 2.2 g / cm 3 without perturbation Indicated.
  • FIG. 2A shows the recording of the horizontal particle velocity field Vx
  • FIG. 2B shows the recording of the vertical particle velocity field Vz, respectively.
  • the P and S waves in the horizontal particle velocity field change their polarity at the source position when the vertical impulse source is used. On the contrary, as shown in FIG. You can see that does not change.
  • the present inventors have shown an extended order propagator calculated as an EGS propagator according to an embodiment of the present invention implemented as described above in a medium in which microvariation exists.
  • FIG. 3 is a view showing a wave field propagated to an EGS propagator according to an embodiment of the present invention in a medium in which microvariation exists.
  • FIG. 3A is a wave field of P wave when given a vertical point sound source on a model given a constant microvariate
  • FIG. 3B is a wave of S wave when given a vertical point sound source on a model given a constant microvariate Chapters are shown respectively.
  • each of the P and S wave fields represents a wave field divided into P and S waves by [M 0 ] -1 of the above algorithm, and the P wave speed and the S wave speed in this model are 2100 m / s, respectively. and 700 m / s, and the density is 2.2 g / cm 3, and the velocity and density of the micro-variance (perturbation) approximation (approximation) background rate of the S wave of the P wave and the medium in order to determine the degree of the waveform when the given is 2/3 of the actual speed was given.
  • the elastic split-step represents zero-order expansion
  • EGSP1 represents a first order expansion
  • EGSP2 represents a second order extended propagator proposed by the present invention.
  • FIG. 4 is a P-wave with a wave length and a mode separation operator for each propagated component propagated to an EGS propagator according to an embodiment of the present invention in a simple two-layer horizontal model. It is a figure which shows the wave field which isolate
  • the present inventors confirmed on the simple horizontal two-layer model whether the wave field is accurately separated into P-wave and S-wave by [M 0 ] -1 , which is a mode separation operator included in the EGS propagator.
  • the P wave speed, S wave speed and density of the upper layer of the used model are 2500 m / s, 1200 m / s and 2.2 g / cm 3 , respectively, and the P wave speed, S wave speed and density of the lower layer are 3500 m / s, 1500 m / s and 2.4 g / cm 3 , the size of the model is 3 ⁇ 2 km and the depth of the upper layer is 0.8 km.
  • FIG. 4 when the horizontal point source is generated at the 1.5 km surface, the wave fields in the vertical particle velocity field and the horizontal particle velocity field are shown in FIG.
  • the wave field divided into P wave and S wave is shown in FIG. 4B, and the picture on the left panel is a snapshot of the wave field 0.5 seconds after the sound source is generated, and the picture on the right panel shows whether the sound source is generated.
  • the wave field after 0.8 second is shown.
  • the mode switching waves generated at the interface that is, SP and PS also exist only in the P and S wave wave fields, that is, the wave fields propagated to the EGS propagator according to the embodiment of the present invention are P and S waves in each propagation step.
  • cross-correlation of a wave field generated by a propagator in a sound source region and a wave field generated by a propagator in a receiver region ie, an input multicomponent wave field used for structural correction
  • imaging is accomplished by convolution, which means that the final structural correction cross section is obtained.
  • the question of how to correlate or convolution is an imaging condition, and in general, if the P wave and the S wave are correctly separated from the acquired multicomponent seismic data, the following equation (18) is used.
  • a structural correction method using an acoustic wave equation can also produce good results.
  • the P wave is generated in the sound source wave field and the S-wave wave field is completely separated from the multicomponent input data in the receiver wave field. Can be obtained.
  • the P wave and the S wave are separated and stored separately in the propagation step by the mode separation operator ([M 0 ] -1 ) in the EGS propagator.
  • the mode separation operator [M 0 ] -1 ) in the EGS propagator.
  • mode separation and imaging are performed in the same frequency-spatial domain, eliminating the need for additional Fourier transformations, and freeing high-quality images from aliasing or noise that can be caused by performing numerical FFT. You can get it.
  • the EGS structure correction algorithm according to the embodiment of the present invention does not use the above Equation (18), but for more stable imaging Schleicher et al.
  • the stabilized division method proposed by 2008 see Ref. 15
  • imaging conditions as shown in Equation 19 below were used.
  • ⁇ shown in the denominator is a value obtained by applying a scaling fraction (0 ⁇ ⁇ 1) to the square of the absolute value of the oscillator wave field.
  • i and j represent the vector wavefield of the sound source and the oscillator (ie, the x-component and z-component in the multi-component seismic data).
  • FIG. 5 is a diagram illustrating a concept of an EGS structure correction algorithm using the EGS propagator illustrated in FIG. 1.
  • FIG. 5 the conceptual diagram of the EGS structure correction algorithm shown in FIG. 5 may be referred to as the case where the thin film propagation described with reference to FIG. 1 is extended to the entire model.
  • FIG. “Indicates wave field propagation in the sound source region, and upward propagation, ie”
  • Backward propagation indicates wave field propagation in the receiver region.
  • each mode separation and mode coupling is performed before and after each thin film propagation, and the P- and S-waves separated from each other are stored in separate arrays and finally used for imaging conditions.
  • Equation 19 in the middle of FIG. 5 represents the correlation of the imaging conditions, and at this time, Equation 19 described above was used as the imaging condition.
  • FIG. 6 is a flowchart schematically illustrating the overall configuration of the EGS structure correction algorithm according to the embodiment of the present invention, which is implemented based on the conceptual diagram illustrated in FIG. 5.
  • the structural correction algorithm is divided into a plurality of components of a seismic wave to be analyzed, establishes a model for a sound source and a receiver, and performs a Fourier transform. Determining a frequency band to be calculated through the step, and using the above EGS radio wave propagation wave propagation from source and wave field in the receiver region for each frequency band Computing a backward propagation from receiver, and integrating the results calculated in the above step through cross correlation, and outputting the image data corrected by the imaging conditions described above. It can be configured to include a series of steps.
  • the step of calculating the wave field propagation in the sound source region is performed by the mode separation operator ([M 0 ] -1 ) in the EGS propagator according to the embodiment of the present invention. separating the wavefield, calculating screen and mode coupling in the fx region, and performing Fourier transform on the spatial variable (x) in the calculated results, then fk Calculating an extrapolated wave field in the region, storing each mode for migration, and reconstructing the sound source wave field by the mode combining operator (M 0 ) in the EGS propagator described above. ) And performing an inverse Fourier transform on the reconstructed wave field.
  • the step of calculating the wave field propagation in the oscillator region likewise, the step of separating the receiver wavefield by the mode separation operator ([M 0 ] -1 ), Computing the screen and mode coupling in the fx region, performing Fourier transform on the spatial variable (x) in the calculated results, and then extrapolating the wave in the fk region Calculating the field, storing each mode for migration, recomposing the oscillator wave field by the mode combining operator M 0 in the EGS propagator, and reconstructed wave field And performing an Inverse Fourier Transform.
  • EGS propagator described above may use the EGS propagator obtained as described with reference to Equations 1 to 17 above.
  • imaging conditions described above may use imaging conditions as described above with reference to Equations 18 to 19.
  • each code is optimized to a Message Passing Interface (MPI) code to enable parallel processing.
  • MPI Message Passing Interface
  • a flowchart for the MPI may be configured as shown in FIG.
  • FIG. 7 is a flowchart schematically illustrating an overall configuration of a process of implementing MPI for parallel processing of an EGS structure correction algorithm according to an embodiment of the present invention illustrated in FIG. 6.
  • the S wave propagates the medium and its polarity is changed at the reflection point of the medium interface, and the conventional elastic structural correction techniques have to separate the S wave polarity reversal phenomenon in the spatial domain.
  • the conventional elastic structural correction techniques have to separate the S wave polarity reversal phenomenon in the spatial domain.
  • the EGS technique according to the embodiment of the present invention can correct the polarity inversion phenomenon in the frequency-wavelength region due to the characteristics of the radio wave calculation.
  • the reflection angle at the reflection point in the waveguide region can be obtained by the following Equation (20).
  • ⁇ represents the reflection angle at the reflection point
  • k h and k z represent the wave number in the distance direction and the wave number in the depth direction, respectively.
  • the sign of the reflection angle at the reflection point can be represented by a sign of horizontal wavenumber of ⁇ 1 ⁇ , so that the polarity of the S wave is in the frequency-wavelength region.
  • k x 0.
  • FIG. 8 is a diagram showing a comparison of structural correction results of applying and not applying the polarity inversion correction method using the imaging conditions shown in Equation 20, respectively.
  • FIG. 8A illustrates a result of performing structural correction through an EGS structure correction technique according to an embodiment of the present invention using synthetic seismic data generated by a single sound source in a simple two-layer model.
  • 8B shows a PP cross section
  • FIG. 8C shows a PS cross section without correcting the polarity inversion phenomenon
  • FIG. 8D shows a PS cross section correcting the polarity inversion phenomenon by using Equation 20 described above.
  • the present inventors applied the above-described EGS structure correction technique to a more realistic and complicated model to confirm whether the image obtained through the actual model and the structure correction fits well, and for this purpose, a benchmark model in the field of seismic exploration In the SEG / EAGE Salt Model (Aminzadeh et al., 1997, see Ref. 16), which is widely known as the SEG / EAGE Salt Model, artificially generates multi-component seismic data through numerical modeling, and uses the same according to an embodiment of the present invention.
  • the structural correction section was obtained by applying the EGS structural correction technique.
  • FIG. 9 is a diagram showing the P wave velocity structure of the two-dimensional cross-sectional model generated as described above.
  • sound sources were installed at intervals of 40 m in the 760 to 12760 m area on the surface of the model, and the receiver was installed at intervals of 20 m in the 60 to 13460 m area to acquire the synthesized seismic data. Is 12.5 Hz.
  • FIG. 10 is a diagram illustrating images of final PP and PS cross sections derived by applying the EGS structure correction technique according to an exemplary embodiment of the present invention.
  • FIG. 10A is a result of the PP structure correction, which shows that the image of the salt dome is represented very accurately compared to the velocity image of FIG. 9, and the peripheral reflection boundary is also imaged at the correct position. You can check it.
  • FIG. 10B is a PS result of performing structure correction without correcting the S-wave polarity inversion phenomenon
  • FIG. 10C is a PS result of using the polarity inversion correction module included in the EGS structure correction algorithm according to an embodiment of the present invention.
  • the image quality is deteriorated when the polarity inversion is not corrected, and in particular, in the case of the Salt upper event, it is shown in the form of a dot rather than a smooth curve. Long horizontal events do not appear well in results without polarity inversion compensation.
  • the EGS propagator and the EGS structure correction method using the same according to the embodiment of the present invention are more accurate and realistic structure correction and imaging than the methods of the prior art to improve the quality of the structure correction image It can be seen that.
  • the vertical slow term is provided by providing the pre-polymerization EGS structure correction method for the multi-component seismic data, which is configured to realize more accurate wave propagation even in horizontally varying speed models. It is possible to solve the problem of the prior art EGS propagator, which had a limitation of approximating to the second term.
  • a P wave and a shot wave is directly used as a structure correction input without having to separate the input multicomponent wave field into P and S waves, including a PS separation module implemented inside the EGS propagator.
  • EGS structure correction method is provided for pre-polymerization of ESA multi-component data configured to generate S-wave image cross-sections. Computation and structure are complicated by separating wave field into P wave and S wave as input before structure correction. It is possible to solve the problems of the conventional structure correction techniques, which had disadvantages.
  • the present invention before the S-wave imaging, by adding a module for correcting the polarity switching in the frequency-wavelength region to improve the quality of the S-wave structure correction images
  • polarity is reversed at the reflection point in the case of seismic S-wave collection, so if it is used for structure correction as it is, the continuity of the event is reduced and the quality of the final cross section is reduced. It is possible to solve the problems of the structural correction techniques of the prior art that had a problem of deterioration.

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geophysics (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • Theoretical Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Software Systems (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Geometry (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • Computer Graphics (AREA)
  • Computational Mathematics (AREA)
  • Algebra (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Multimedia (AREA)
  • Pure & Applied Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Electromagnetism (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Radar Systems Or Details Thereof (AREA)
  • Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)

Abstract

본 발명은 수평적 속도 변화가 존재하는 매질에서 파의 전파를 빠르게 계산할 수 있는 기존의 SGS(Scalar Generalized-Screen) 기법을 탄성 파동방정식(elastic wave equation)으로 확장하는 것에 의해, 지하 매질의 경계면들을 전파하면서 P파와 S파 상호간의 모드 전환을 거치는 탄성파의 거동을 효율적으로 표현할 수 있는 EGS(Elastic Generalized-Screen) 전파자(wave propagator)를 이용한 단방향 파동방정식(one-way wave equation) 중합 전 심도 구조보정(prestack depth migration) 방법에 관한 것으로, 본 발명에 따르면, 전파자의 수직 느리기항(vertical slowness term)의 테일러 급수전개를 2차까지 확장하여 복잡한 구조의 매질에서 보다 높은 정확도의 파동장을 계산할 수 있으며, 모드 분리 연산자를 전파자에 포함시킴으로써 다성분 자료를 P파와 S파로 분리할 필요 없이 음원 모음(shot gather)을 바로 구조보정 입력으로 사용하여 P파 및 S파 영상 단면을 생성할 수 있는 데 더하여, S파 영상화 전에 파수-주파수 영역에서 극성 전환에 대한 보정을 행하여 S파 구조보정 영상의 품질을 향상시킬 수 있도록 구성되는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법이 제공된다.

Description

탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법
본 발명은 탄성파 구조보정(migration) 방법에 관한 것으로, 더 상세하게는, 수평적 속도 변화가 존재하는 매질에서 파의 전파를 빠르게 계산할 수 있는 기존의 스칼라 일반화된 막(Scalar Generalized-Screen ; SGS) 기법을 탄성 파동방정식(elastic wave equation)으로 확장하는 것에 의해, 지하 매질의 경계면들을 전파하면서 P파와 S파 상호간의 모드 전환을 거치는 탄성파의 거동을 효율적으로 표현할 수 있는 탄성 일반화된 막(Elastic Generalized-Screen ; EGS) 전파자(wave propagator)를 이용한 단방향 파동방정식(one-way wave equation) 중합 전 심도 구조보정(prestack depth migration) 방법에 관한 것이다.
또한, 본 발명은, 전파자(wave propagator)의 수직 느리기항(vertical slowness term)의 테일러 급수전개를 2차까지 확장하여 복잡한 구조의 매질에서 보다 높은 정확도의 파동장을 계산할 수 있도록 구성되는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법에 관한 것이다.
아울러, 본 발명은, 다성분 자료의 구조보정시 입력자료를 P파와 S파 파동장으로 먼저 분리한 후 사용하는 종래의 기법들과 달리, 모드 분리 연산자를 전파자에 포함시킴으로써 다성분 자료를 P파와 S파로 분리할 필요 없이 음원 모음(shot gather)을 바로 구조보정 입력으로 사용하여 P파 및 S파 영상 단면을 생성할 수 있도록 구성되는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법에 관한 것이다.
더욱이, 본 발명은, 탄성파 S파 음원 모음(shot gather)의 경우 반사지점에서 극성이 역전되는 현상이 있음으로 인해 이를 그대로 구조보정에 사용하게 되면 이벤트의 연속성이 떨어지게 되어 최종 단면의 품질이 저하되는 종래기술의 문제점을 해결하기 위해, S파 영상화 전에 파수-주파수 영역에서 극성 전환에 대한 보정을 수행하는 모듈을 추가하여 S파 구조보정 영상의 품질을 향상시킬 수 있도록 구성되는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법에 관한 것이다.
일반적으로, 탄성파 구조보정(migration)이란, 중합 단면(stack section) 상의 모든 일차 반사 이벤트들을 제 위치로 옮겨 주고, 회절(diffraction) 곡선을 중합(collapse) 하여 공간 분해능을 높이고 지하구조의 영상을 얻도록 해주는 과정을 의미한다.
여기서, 보정 위치는 반사파의 주행시간(travel time)과 속도(velocity)에 의해 결정되며, 이러한 구조보정 방법에는, 크게 나누어, 파동 방정식의 적분 해에 기초한 키르히호프(Kirchhoff) 구조보정, 유한 차분(finite difference)법을 이용한 구조보정 및 주파수와 파수를 이용한 주파수-파수(F-K) 구조보정으로 나누어진다.
더 상세하게는, 먼저, 키르히호프(Kirchhoff) 구조보정법은, 파선 이론을 기초로 하여 중합 단면상의 각 점에서 그 점을 정점으로 하는 적절한 회절 곡선(diffraction hyperbola)을 계산하고, 이 회절 곡선상에 놓이는 모든 샘플값을 합하여 구조보정 단면상에 대응하는 점의 진폭으로 취하는 방법으로, 계산의 속도가 빠르고 비교적 정확한 결과를 도출하여 산업계에서는 많이 사용되는 방법이나, 복잡한 구조의 경우 파선 이론이 제대로 적용되지 못하여 실패할 수도 있는 단점이 있다.
또한, 유한 차분을 이용한 구조보정 기법에는, 일반적으로, 양방향 파동 방정식(full two-way wave equation)을 사용한 역시간 구조보정(reverse time migration) 방법이 가장 널리 사용되는데, 이는 복잡한 구조에서도 비교적 정확한 결과를 도출하는 장점이 있지만 계산시간이 너무 많이 걸리는 단점이 있다.
아울러, 단방향 파동 방정식(one-way wave equation)을 이용한 구조보정 기법은, 상기한 바와 같이 복잡한 구조에서 영상화에 실패하는 경우가 있는 키르히호프(Kirchhoff) 방식의 단점과, 계산시간이 많이 걸리는 역시간 구조보정의 단점을 보완하여 제시된 방법으로, 양방향 파동 방정식을 근사시켜 한쪽 방향으로만 전파하는 상향파(upgoing wave)만을 고려하도록 고안되어 시간 영역에서 유한 차분으로 풀거나 주파수-파수(F-K) 영역으로 옮겨서 파를 전파시키는 방법이다.
이때, 파를 전파시키는 전파자(propagator)로서 단방향 파동 방정식을 사용하고, 음원 영역에서 발생한 파동장과 수진기 영역에서 발생한 파동장을 영상화 조건(imaging condition)에 적용하여 최종 구조보정 단면을 얻는 것이 단방향 파동 방정식 구조보정 기법이다.
상기한 바와 같이, 현재까지 대부분의 탄성파 탐사자료의 영상화(imaging)는 스칼라 파동 방정식(scalar wave equation)을 기반으로 수행되어 왔다.
여기서, 해양에서 스트리머를 이용하여 취득한 탄성파 자료나 육상에서 단성분 지오폰(geophone)으로 취득한 자료의 경우는 이러한 파동 방정식으로 처리가 가능하였으나, 지구의 지하 매질은 엄밀히 말하면 물과 같은 음향(acoustic) 매질이 아닌 복잡하고 비균질(heterogeneous) 특성을 가진 탄성 매질(elastic media)이다.
즉, 종래에는, 이러한 탄성 매질에서 취득한 자료를 스칼라 파동 방정식(scalar wave equation) 또는 음향 파동 방정식(acoustic wave equation)으로 가정하여 처리하는 방법을 사용해 왔으나, 이는, 탄성파(elastic wave)를 음향파(acoustic wave)로 간주하는 방법이므로 엄밀하게 말하면 정확한 방법이 아니다.
따라서 탄성 매질을 진행하여 돌아온 탄성파는 그 거동을 정확히 묘사할 수 있도록 매질의 수평 및 수직 변위(또는 속도, 가속도)를 탐지할 수 있는 3성분(3-component)의 지오폰(geophone)으로 탐지하여야 하며, 그 자료는 탄성 파동 방정식(elastic wave equation)을 이용하여 처리해야 한다.
또한, 최근에는, 탐사장비와 컴퓨팅 환경의 발달로 인해 석유탐사 시장 등에서 다성분(multi-component) 탄성파 탐사자료를 많이 취득하고 있으며, 여기에 탄성 파동 방정식(elastic wave equation)을 이용한 구조보정을 많이 실시하고 있다.
즉, 탄성 파동 방정식을 이용한 구조보정 기법은, 상기한 스칼라 파동 방정식을 이용한 구조보정 기법과 그 방식은 비슷하나, 입력 파동장으로 다성분 탄성파 자료가 사용되며, 파동 방정식 또한 음향 파동 방정식이 아닌 탄성 파동 방정식을 사용하여 탄성파가 매질을 진행하면서 발생하는 모드의 전환(P파에서 S파 또는 S파에서 P파로 모드가 전환되는 현상)과 각 모드의 파(즉, P파와 S파)의 감쇠를 효율적으로 표현할 수 있다.
더 상세하게는, 종래, 예를 들면, Hokstad(2000)나 Sun and McMechan(2011) 등은 Kirchhoff 구조보정 기법을 탄성(elastic) 으로 확장하여 다성분 탐사 자료에 적용하였으며, Chang and McMechan(1994)과 Yan and Sava(2008)는 탄성 파동 방정식을 이용하여 역시간 구조보정 기법을 개발하였다(참고문헌 1 내지 참고문헌 4 참조).
그러나 스칼라 파동방정식 기법들과 마찬가지로, 탄성(Elastic) Kirchhoff 기법은 여전히 복잡한 구조에서는 영상화에 실패를 하는 경우가 발생하고, 탄성(Elastic) 역시간 구조보정 기법은 계산 비용이 많이 드는 단점이 존재한다.
이에, 이러한 종래기술의 단점들을 보완하는 방법으로, 탄성(Elastic) 단방향 파동 방정식이 제시되었으며, 즉, Landers and Claerbout(1972)는 처음으로 위상-막(phase-screen) 기법을 음향(acoustic) 파동방정식에서 탄성(elastic) 파동방정식으로 확장하였고, Fisk and McCartor(1991)는 P파와 S파의 모드 전환을 단방향 파동 방정식으로 구현하였으며, Xie and Wu(2005)는 위상-막(phase-screen) 기법을 확장하여 구조보정 모듈을 개발하였으나, 이들은 입력 파동장으로 모드가 분리된 P파와 S파를 각각 사용하였고 그들의 구조보정 모듈에는 모드 전환이 고려되지 않았다(참고문헌 5 내지 참고문헌 7 참조).
또한, Le Rousseau and De Hoop(2003)는 F-K 구조보정 기법인 스플릿-스텝(split-step) 기법과 위상-막(phase-screen) 기법을 탄성파(elastic wave)에 일반화(generalization)하여 탄성 일반화된 막(elastic generalized-screen ; EGS) 전파자(propagator)를 제시한 바 있다(참고문헌 8 참조).
그러나 이들 방법은, 수직 느리기 항(vertical slowness operator)을 1차항 까지만 근사하였고, 전파자 개발 단계에서 멈추어 더 이상 구조보정 기법으로는 개발되지 않은 점에서 한계가 있는 것이었다.
따라서 상기한 바와 같은 종래기술의 문제점을 해결하기 위하여는, 수직 느리기 항이 1차항까지 밖에 계산되지 않았던 종래의 EGS 전파자를 2차항까지 확장하여 복잡한 모델이나 수평적으로 속도 변화가 심한 모델에서도 보다 정확한 파의 전파를 구현하는 것에 의해 S파의 극성 반전 현상을 보정할 수 있도록 구성되어 EGS 전파자의 성능을 더욱 향상시킬 수 있는 새로운 구조보정 알고리즘을 제공하는 것이 바람직하나, 아직까지 그러한 요구를 모두 만족시키는 장치나 방법은 제시되지 못하고 있는 실정이다.
[선행기술문헌]
1. 한국 등록특허공보 제10-1413751호 (2014.06.24.)
2. 한국 등록특허공보 제10-1219746호 (2013.01.02.)
3. 한국 등록특허공보 제10-1347969호 (2013.12.27.)
4. 한국 등록특허공보 제10-1281803호 (2013.06.27.)
[참고문헌]
1. Hokstad, K., 2000, Multicomponent Kirchhoff migration, Geophysics, 65, 861-873.
2. Sun R., and G. A. McMechan, 2011, Prestack 2D parsimonious Kirchhoff depth migration of elastic seismic data, Geophysics, 76, s157-s164.
3. Chang, W. F., and McMechan, G. A., 1994, 3-D elastic prestack reverse-time depth migration, Geophysics, 59, 579-609.
4. Yan J., and Sava P., 2008, Isotropic angle-domain elastic reverse-time migration, Geophysics, 77, 229-239.
5. Landers, T., and Claerbout, J. F., 1972, Numerical calculation of elastic waves in laterally inhomogeneous media, J. Geophys. Res., 77, 1476-1482.
6. Fisk, M. D., and McCartor, G. D., 1991, The phase screen method for vector elastic waves, J. Geophys. Res., 96, 5985-6010.
7. Xie X. B., and Wu, R. S., 2005, Multicomponent prestack depth migration using the elastic screen method, Geophysics, 70, 30-37.
8. Le Rousseau, J. H. and De Hoop M. V., 2003, Generalized-screen approximation and algorithm for the scattering of elastic waves, Q. JI Mech. Appl. Math., 56, 1-33.
9. Le Rousseau, J. H., 2001, Microlocal analysis of wave-equation imaging and generalized-screen propagators, Ph. D. Thesis, CWP, CSM.
10. Stoffa, P. L.,Fokkema, J. T., De Luna F. R. M., and Kessinger, W. P., 1990, Split-step Fourier Migration, Geophysics, 55, 410-421.
11. Wu, R. S., and L. J. Huang, 1992, Scattered field calculation in heterogeneous media using phase-screen propagator, Expanded Abstracts of SEG 1992 Annual Meeting, 1289-1292.
12. De Hoop, M. V., 1996, Generalization of the Bremmer coupling series, J. Math. Phys., 37, 3246-3282.
13. De Hoop, M. V., and De Hoop, A. T., 1994, Elastic wave up/down decomposition in inhomogeneous and anisotropic media: an operator approach and its approximations, Wave Motion, 20, 57-82.
14. Sava P. C. and Fomel S., 2003, Angle-domain common-image gathers by wavefield continuation method, Geophysics, 68, 1065-1074.
15. Schleicher J., Costa J. C, and Novais A., 2008, A comparison of imaging condition for wave-equation shot-profile migration, Geophysics, 73, S219-S227.
16. Aminzadeh, F., J. Brac, and T. Kunz, 1997, SEG/EAGE 3-D salt and overthrust models, in SEG/EAGE 3-D Modeling Series, No.1, SEG.
본 발명은 상기한 바와 같은 종래기술의 문제점을 해결하고자 하는 것으로, 따라서 본 발명의 목적은, 수직 느리기 항(vertical slowness operator)이 1차항 까지만 근사되는 한계가 있었던 종래기술의 EGS 전파자의 문제점을 해결하기 위해, 수직 느리기 항을 2차항까지 확장하여 복잡한 모델이나 수평적으로 속도 변화가 심한 모델에서도 보다 정확한 파의 전파를 구현하여 EGS 전파자의 성능을 더욱 향상시킬 수 있도록 구성되는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법을 제공하고자 하는 것이다.
또한, 본 발명의 다른 목적은, 구조보정 전에 파동장을 P파와 S파로 분리하여 입력으로 사용함으로 인해 계산 및 구조가 복잡해지는 단점이 있었던 종래기술의 구조보정 기법들의 문제점을 해결하기 위해, EGS 전파자 내부에 구현된 P-S분리 모듈을 포함하여 입력 다성분 파동장을 P파와 S파로 분리할 필요 없이 음원 모음(shot gather)을 바로 구조보정 입력으로 사용하여 P파 및 S파 영상 단면을 생성할 수 있도록 구성되는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법을 제공하고자 하는 것이다.
아울러, 본 발명의 또 다른 목적은, 탄성파 S파 음원 모음(shot gather)의 경우 반사지점에서 극성이 역전되는 현상이 있음으로 인해 이를 그대로 구조보정에 사용하게 되면 이벤트의 연속성이 떨어지게 되어 최종 단면의 품질이 저하되는 문제가 있었던 종래기술의 구조보정 기법들의 문제점을 해결하기 위해, S파 영상화 전에 주파수-파수 영역에서 극성 전환에 대한 보정을 수행하는 모듈을 추가하여 S파 구조보정 영상의 품질을 향상시킬 수 있도록 구성되는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법을 제공하고자 하는 것이다.
상기한 바와 같은 목적을 달성하기 위해, 본 발명에 따르면, 지하 매질의 경계면들을 전파하면서 P파와 S파 상호간의 모드 전환을 거치는 탄성파의 거동을 표현하기 위한 EGS(Elastic Generalized-Screen) 전파자(wave propagator)의 수직 느리기항(vertical slowness term)을 2차항까지 확장하여 다성분 자료의 구조보정시 입력자료를 P파와 S파 파동장으로 분리할 필요 없이 음원 모음(shot gather)을 바로 구조보정 입력으로 사용하여 P파 및 S파 영상 단면을 생성할 수 있도록 구성되는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법에 있어서, 분석하고자 하는 탄성파 다성분 자료를 입력받아 음원(source)과 수진기(receivers)에 대한 모델을 수립하고, 푸리에 변환(Fourier transform)을 통하여 계산해야 할 주파수 대역(frequency band)을 결정하는 전처리 단계; 상기 EGS 전파자를 이용하여 각각의 주파수 대역에 대하여 음원 영역에서의 파동장 전파(forward propagation from source)를 산출하는 음원영역 처리단계; 상기 EGS 전파자를 이용하여 각각의 주파수 대역에 대하여 수진기 영역에서의 파동장 전파(backward propagation from receiver)를 산출하는 수진기영역 처리단계; 상기 음원영역 처리단계 및 상기 수진기영역 처리단계에서 각각 계산된 결과를 상호상관(Cross Correlation)을 통해 통합하고 영상화 조건(Imaging condition)에 의해 구조보정을 행하는 구조보정단계; 및 상기 구조보정단계에서 구조보정된 영상 데이터를 출력하는 출력단계를 포함하여 구성되는 것을 특징으로 하는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법이 제공된다.
여기서, 상기 EGS 전파자는, 이하의 수학식에 나타낸 EGS 전파자를 이용하도록 구성되는 것을 특징으로 한다.
Figure PCTKR2015010951-appb-I000001
(여기서, xμ(μ = 1, 2)과 x3는 각각 수평과 수직 좌표를 나타내고, s = -iω(ω는 각주파수) 이며, αv = kv/iω(kv는 수평성분 파수, v = 1, 2)이고,
Figure PCTKR2015010951-appb-I000002
와
Figure PCTKR2015010951-appb-I000003
는 각각 Fourier 변환과 그의 역변환을 나타내며, M0는 고유치벡터 (eigenvector)로 이루어진 대각화 행렬로서 분리된 P파와 S파를 결합하는(coupling) 연산자로 작용하고, [M0]-1은 M0의 역행렬로서 P파와 S파를 분리하는 연산자로 작용하며, λ는 항(term)의 수를 나타냄.)
또한, 상기 음원영역 처리단계는, 상기 EGS 전파자 내의 모드 분리 연산자(operator)([M0]-1)에 의해 음원 파동장(source wavefield)을 분리하는 단계; 주파수-공간(f-x) 영역에서 스크린(screen) 및 모드 결합(mode coupling)을 계산하는 단계; 계산된 결과에 공간 변수(spatial variable, x)에 대하여 푸리에 변환을 수행하는 단계; 주파수-파수(f-k) 영역에서 외삽된(extrapolated) 파동장을 계산하는 단계; 구조보정(migration)을 위해 각각의 모드를 저장하고, 상기 EGS 전파자 내의 모드결합 연산자(M0)에 의해 상기 음원 파동장을 재구성(recomposition)하는 단계; 및 재구성된 상기 음원 파동장에 역푸리에 변환을 수행하는 단계를 포함하여 구성되는 것을 특징으로 한다.
아울러, 상기 수진기영역 처리단계는, 상기 EGS 전파자 내의 모드 분리 연산자(operator)([M0]-1)에 의해 수진기 파동장(receiver wavefield)을 분리하는 단계; 주파수-공간(f-x) 영역에서 스크린(screen) 및 모드 결합(mode coupling)을 계산하는 단계; 계산된 결과에 공간 변수(spatial variable, x)에 대하여 푸리에 변환을 수행하는 단계; 주파수-파수(f-k) 영역에서 외삽된(extrapolated) 파동장을 계산하는 단계; 구조보정(migration)을 위해 각각의 모드를 저장하고, 상기 EGS 전파자 내의 모드결합 연산자(M0)에 의해 상기 수진기 파동장을 재구성(recomposition)하는 단계; 및 재구성된 상기 수진기 파동장에 역푸리에 변환을 수행하는 단계를 포함하여 구성되는 것을 특징으로 한다.
여기서, 상기 음원영역 처리단계 및 상기 수진기영역 처리단계는, 복수의 프로세서를 이용하여, 각각의 주파수 대역별로 처리를 할당하여 병렬처리(parallel processing)가 가능하도록 구성됨으로써, 전체적인 처리시간을 단축할 수 있도록 구성되는 것을 특징으로 한다.
더욱이, 상기 구조보정단계는, 상기 영상화 조건으로서, 이하의 수학식에 나타낸 영상화 조건을 이용하도록 구성되는 것을 특징으로 한다.
Figure PCTKR2015010951-appb-I000004
(여기서, Iij는 최종 영상화 이미지,
Figure PCTKR2015010951-appb-I000005
는 푸리에 영역에서의 스칼라(scalar) 음원 파동장(source wavefield),
Figure PCTKR2015010951-appb-I000006
은 푸리에 영역에서의 스칼라 수진기 파동장(receiver wavefield)을 나타내고, ε는
Figure PCTKR2015010951-appb-I000007
로 표현할 수 있으며, i 와 j 는 음원과 수진기의 벡터 파동장(vector wavefield)으로, 다성분 탄성파 자료에서 수평성분(x-component)과 수직성분(z-component)을 각각 의미함.)
또한, 상기 방법은, 이하의 수학식을 이용하여 주파수-파수 영역에서 반사 지점에서의 반사각을 구하는 것에 의해, 매질 경계면의 반사점에서 극성이 변화하게 되는 S파 극성 역전 현상을 보정하는 극성보정단계를 더 포함하여 구성되는 것을 특징으로 한다.
Figure PCTKR2015010951-appb-I000008
(여기서, γ는 반사지점에서의 반사각을 나타내며, kh와 kz는 각각 거리 방향의 파수와 심도 방향의 파수를 나타냄.)
아울러, 본 발명에 따르면, 상기에 기재된 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법을 이용하여, 다성분 자료를 P파와 S파로 분리할 필요 없이 음원 모음(shot gather)을 바로 구조보정 입력으로 사용하여 P파 및 S파 영상 단면을 생성할 수 있도록 구성되는 동시에, S파 영상화 전에 파수-주파수 영역에서 극성 전환에 대한 보정을 수행하여 S파 구조보정 영상의 품질을 향상시킬 수 있도록 구성되는 것을 특징으로 하는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 시스템이 제공된다.
상기한 바와 같이, 본 발명에 따르면, 수직 느리기 항(vertical slowness operator)을 2차항까지 확장하여 복잡한 모델이나 수평적으로 속도 변화가 심한 모델에서도 보다 정확한 파의 전파를 구현하여 EGS 전파자의 성능을 더욱 향상시킬 수 있도록 구성되는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법이 제공됨으로써, 수직 느리기 항이 1차항 까지만 근사되는 한계가 있었던 종래기술의 EGS 전파자의 문제점을 해결할 수 있다.
또한, 본 발명에 따르면, EGS 전파자 내부에 구현된 P-S분리 모듈을 포함하여 입력 다성분 파동장을 P파와 S파로 분리할 필요 없이 음원 모음(shot gather)을 바로 구조보정 입력으로 사용하여 P파 및 S파 영상 단면을 생성할 수 있도록 구성되는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법이 제공됨으로써, 구조보정 전에 파동장을 P파와 S파로 분리하여 입력으로 사용함으로 인해 계산 및 구조가 복잡해지는 단점이 있었던 종래기술의 구조보정 기법들의 문제점을 해결할 수 있다.
아울러, 본 발명에 따르면, S파 영상화 전에 주파수-파수 영역에서 극성 전환에 대한 보정을 수행하는 모듈을 추가하여 S파 구조보정 영상의 품질을 향상시킬 수 있도록 구성되는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법을 제공됨으로써, 탄성파 S파 음원 모음(shot gather)의 경우 반사지점에서 극성이 역전되는 현상이 있음으로 인해 이를 그대로 구조보정에 사용하게 되면 이벤트의 연속성이 떨어지게 되어 최종 단면의 품질이 저하되는 문제가 있었던 종래기술의 구조보정 기법들의 문제점을 해결할 수 있다.
도 1은 EGS 전파자의 개념을 나타내는 도면이다.
도 2는 균질(zero-pertubation) 매질에서 본 발명의 실시예에 따른 EGS 전파자로 전파한 파동장을 나타내는 도면이다.
도 3은 미소변량이 존재하는 매질에서 본 발명의 실시예에 따른 EGS 전파자로 전파한 파동장을 나타내는 도면이다.
도 4는 단순 2층 수평모델에서 본 발명의 실시예에 따른 EGS 전파자로 전파한 전파시킨 각 성분(component)별 파동장과 모드 분리 연산자(operator)로 P파와 S파를 분리한 파동장을 각각 나타내는 도면이다.
도 5는 도 1에 나타낸 EGS 전파자를 이용한 EGS 구조보정 알고리즘의 개념을 나타내는 도면이다.
도 6은 도 5에 나타낸 바와 같은 개념도를 바탕으로 구현된 본 발명의 실시예에 따른 EGS 구조보정 알고리즘의 전체적인 구성을 개략적으로 나타내는 플로차트이다.
도 7은 도 6에 나타낸 본 발명의 실시예에 따른 EGS 구조보정 알고리즘의 병렬처리를 위한 MPI를 구현하는 과정의 전체적인 구성을 개략적으로 나타내는 플로차트이다.
도 8은 본 발명의 실시예에 따른 EGS 구조보정 알고리즘의 영상화 조건을 이용하여 극성 역전의 보정 방법을 적용한 경우와 적용하지 않은 경우의 구조보정 결과를 각각 비교하여 나타내는 도면이다.
도 9는 본 발명의 실시예에 따른 EGS 구조보정 알고리즘의 검증을 위해 전 생성된 2차원 단면 모델의 P파 속도 구조를 나타내는 도면이다.
도 10은 본 발명의 실시예에 따른 EGS 구조보정 기법을 적용하여 도출한 최종 PP 및 PS 단면의 영상을 나타내는 도면이다.
이하, 첨부된 도면을 참조하여, 본 발명에 따른 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법의 구체적인 실시예에 대하여 설명한다.
여기서, 이하에 설명하는 내용은 본 발명을 실시하기 위한 하나의 실시예일 뿐이며, 본 발명은 이하에 설명하는 실시예의 내용으로만 한정되는 것은 아니라는 사실에 유념해야 한다.
또한, 이하의 본 발명의 실시예에 대한 설명에 있어서, 종래기술의 내용과 동일 또는 유사하거나 당업자의 수준에서 용이하게 이해하고 실시할 수 있다고 판단되는 부분에 대하여는, 설명을 간략히 하기 위해 그 상세한 설명을 생략하였음에 유념해야 한다.
즉, 본 발명은, 후술하는 바와 같이, 수직 느리기 항(vertical slowness operator)이 1차항 까지만 근사되는 한계가 있었던 종래기술의 EGS 전파자의 문제점을 해결하기 위해, 수직 느리기 항을 2차항까지 확장하여 복잡한 모델이나 수평적으로 속도 변화가 심한 모델에서도 보다 정확한 파의 전파를 구현하여 EGS 전파자의 성능을 더욱 향상시킬 수 있도록 구성되는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법에 관한 것이다.
또한, 본 발명은, 후술하는 바와 같이, 구조보정 전에 파동장을 P파와 S파로 분리하여 입력으로 사용함으로 인해 계산 및 구조가 복잡해지는 단점이 있었던 종래기술의 구조보정 기법들의 문제점을 해결하기 위해, EGS 전파자 내부에 구현된 P-S분리 모듈을 포함하여 입력 다성분 파동장을 P파와 S파로 분리할 필요 없이 음원 모음(shot gather)을 바로 구조보정 입력으로 사용하여 P파 및 S파 영상 단면을 생성할 수 있도록 구성되는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법에 관한 것이다.
아울러, 본 발명은, 후술하는 바와 같이, 탄성파 S파 음원 모음(shot gather)의 경우 반사지점에서 극성이 역전되는 현상이 있음으로 인해 이를 그대로 구조보정에 사용하게 되면 이벤트의 연속성이 떨어지게 되어 최종 단면의 품질이 저하되는 문제가 있었던 종래기술의 구조보정 기법들의 문제점을 해결하기 위해, S파 영상화 전에 주파수-파수 영역에서 극성 전환에 대한 보정을 수행하는 모듈을 추가하여 S파 구조보정 영상의 품질을 향상시킬 수 있도록 구성되는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법에 관한 것이다.
즉, 일반적으로, 전파자(Progagator)라 함은, 지구모델을 수치적으로 설정한 후 그 모델 상에서 탄성파가 진행하는 것을 수치적으로 모사(simulation) 할 때 사용되는 엔진과 같은 연산자(operator) 또는 함수로서, 따라서 전파자가 잘 구현될수록 수치 모델링에 있어 실제 탄성파가 지구 내부를 전파하는 것을 잘 묘사할 수 있다고 볼 수 있다.
여기서, 종래, 예를 들면, Le Rousseau and De Hoop. 2003(참고문헌 8 참조)에서 제시된 전파자는 수직 느리기 항(실제 파가 하부로 전파할 때 사용되는 함수)이 1차 오더(order) 까지만 계산된 것으로, 쉽게 말하면 부정확한 전파자에 해당하고, 즉, 상기한 종래기술 문헌에 제시된 전파자는 단지 수학적으로 전파자만 도출하는 것에 그치는 것일 뿐, 구조보정 알고리즘으로 구현된 바는 없었다.
이에, 본 발명자들은, 1차 오더까지 계산된 종래의 전파자를 개선하여 구조보정 알고리즘으로 구현하는 동시에, 보다 정확한 파의 전파를 위해 EGS 전파자를 2차 오더로 확장하여 수치적으로 전파자의 성능을 향상시켰으며, 이와 같이 하여 향상된 전파자를 이용하여 가지고 구조보정 알고리즘을 구현하였고, 이때, 전파자에서 P파와 S파를 분리하여 구조보정에 사용할 수 있도록 하는 동시에, S파의 극성 역전 현상을 보정하는 모듈을 추가하여 EGS 전파자의 성능을 더욱 향상시킬 수 있도록 하였다.
계속해서, 도면을 참조하여, 본 발명에 따른 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법의 구체적인 구성에 대하여 설명한다.
여기서, 본 발명에 따른 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법의 구체적인 구성에 대하여 설명하기 전에, 먼저, 탄성파 얇은-막(Elastic Thin-slab) 전파자(Propagator)에 대하여 설명하면, 종래, Le Rousseau, 2001(참고문헌 9 참조)는 푸리에 영역 스플릿 스텝(split-step Fourier)(Stoffa et al., 1990, 참고문헌 10 참조) 기법과 위상-막(phase-screen method)(Wu and Huang, 1992, 참고문헌 11 참조) 기법을 수평적으로 속도 변화가 심한 매질에서 보다 정확하게 전파할 수 있도록 일반화하고 이 방법을 일반화된 막(generalized-screen) 전파자라 명명하였으며, Le Rousseau and De Hoop, 2003(참고문헌 8 참조)는 이를 탄성파(elastic wave)로 확장하였다.
더 상세하게는, 만일 매질의 속성(properties)이 두 수직 심도 [x'3, x3] 사이에서 매우 작게 변하고 △x3( = x3 - x'3) 가 충분히 작다면(thin slab), 아주 얇은 막(thin slab) 사이에서의 상향으로 전파하는 상향(upgoing)파와 하향으로 전파하는 하향(downgoing)파를 나타내는 탄성(elastic) 얇은 막(thin-slab) 전파자는 이하의 [수학식 1]과 같이 의사 미분 연산자(pseudo-differential operator,ΨDO)의 형태로 표현된다.
[수학식 1]
Figure PCTKR2015010951-appb-I000009
여기서, xμ(μ = 1, 2)과 x3는 각각 수평과 수직 좌표를 나타내고, s = -iω(ω는 각주파수) 이며, αv = kv/iω(kv는 수평성분 파수, v = 1, 2)이다.
의사 미분 연산자로 표현된 [수학식 1]은 주파수-파수 영역에서 매질의 수직느리기항(γ)을 이용하여 위상을 이동시킴으로써 심도 방향으로 파동장을 외삽 (extrapolation) 해주는 연산자로 공간 변수(spatial variable, x)와 Fourier 변수 (αv)를 하나의 항으로 표현하게 된다.
이때,
Figure PCTKR2015010951-appb-I000010
와
Figure PCTKR2015010951-appb-I000011
는 각각 Fourier 변환과 그의 역 변환을 나타내는 커널이다.
또한, 수평적으로 속도변화가 있는 매질에서의 수직 느리기 항은 이하의 [수학식 2]와 같이 배경 매질에서의 위상 이동만을 고려한 수직 느리기항
Figure PCTKR2015010951-appb-I000012
과 수평적 속도 변화를 고려한 미소 변량(perturbation)항
Figure PCTKR2015010951-appb-I000013
으로 표현된다.
[수학식 2]
Figure PCTKR2015010951-appb-I000014
여기서, 상기한 [수학식 2]에 있어서, 아래첨자 γ는 의사 미분 연산자의 심볼(right symbol)이다.
즉, 이러한 탄성(elastic) 단방향 얇은 막 전파자의 공간영역 변수와 파수영역 변수를 근사에 의해 분리해 내는 것이 수평적 속도 변화를 고려한 탄성 일반화된 막(elastic generalized-screen) 전파자의 핵심이다.
다음으로, 수직 느리기 항의 2차항까지 확장된 일반화된 막(Generalized-screen) 전파자에 대하여 설명한다.
수직 느리기항 연산자(Γ)의 심볼(right symbol), 즉, γr(Le Rousseau and De Hoop, 2003, 참고문헌 8 참조)는 특성화 연산자(characteristic operator, A)에 대하여 이하의 [수학식 3]에 나타낸 바와 같은 관계를 가진다.
[수학식 3]
Figure PCTKR2015010951-appb-I000015
여기서, A의 심볼인 ar을 첫째 항까지만 근사하면(high-frequency approximation)(De Hoop, 1996, 참고문헌 12 참조), 이하의 [수학식 4]에 나타낸 바와 같다.
[수학식 4]
Figure PCTKR2015010951-appb-I000016
여기서, 아래 첨자는 심볼 오더(symbol order)를 나타내고 위첨자는 매질 대비(medium contrast)의 차수를 나타낸다.
아울러, 등방성 매질(isotropic media)에서의 Hooke의 법칙과 De Hoop and De Hoop's(1994) 의 특성 연산자 행렬(characteristic operator matrix)(a)의 심볼 유도식으로부터, 이하의 [수학식 5]와 같이 하여 주 심볼 행렬(principal symbol matrix)(a[2])이 계산될 수 있다.
[수학식 5]
Figure PCTKR2015010951-appb-I000017
여기서, α22, α12, α13 및 α32는 α11, α21, α23 및 α31로부터 α1과 α2를 각각 교환(interchange)함으로써 얻어질 수 있으며, 테일러 급수를 3차항까지 이용하여 ks -1, ks -1/2, ρ1 /2, ρ-1/2를 확장하고(expanding) 이들을 상기한 수학식으로 대체하여 입실론(ε) 차수(order)에 대하여 다시 정리하는(rearrange) 것에 의해, 더욱 고차항의 심볼이 산출될 수 있다.
따라서 얇은 막[x'3, x3] 내에서 배경매질, 즉,
Figure PCTKR2015010951-appb-I000018
의 미소변량(perturbation)은 이하의 [수학식 6]과 같이 표현할 수 있다.
[수학식 6]
Figure PCTKR2015010951-appb-I000019
또한, 주 심볼 행렬(principal symbol matrix)(a[2])에 있는 ks -1, ks -1/2, ρ1 /2 및 ρ-1/2를 테일러 급수로 전개하고 이를 [수학식 6]과 함께 상기한 [수학식 5]에 대입하면, 이하의 [수학식 7] 내지 [수학식 9]에 나타낸 바와 같이 하여 2차항까지 확장된 매질 대비(medium contrast)에 대한 주 심볼(principal symbol) (a)을 얻을 수 있다.
즉, 이하의 [수학식 7]은 상기한 [수학식 5]에서 입실론(ε)에 독립인 항들이며, 따라서 이들은 전파자의 배경항(background term)을 포함한다.
아울러, 이하의 [수학식 7]에 있어서, αγ,22, αγ,12, αγ,13 및 αγ,32는 α1과 α2를 각각 αγ,11, αγ,21, αγ,23 및 αγ,31로 교환함으로써 얻어질 수 있다.
[수학식 7]
Figure PCTKR2015010951-appb-I000020
[수학식 8]
Figure PCTKR2015010951-appb-I000021
[수학식 9]
Figure PCTKR2015010951-appb-I000022
여기서, 상기한 [수학식 8]은 [수학식 5]의 테일러 전개 후 1차 입실론항(first epsilon oredr term)에 의해 재정리된(rearranged) 것이고, [수학식 9]는 [수학식 5]의 테일러 전개 후 2차 입실론항에 의해 재정리된(rearranged) 것이며, 이때, 특정 항들은 루소의 유도식(Le Rousseau's deviation) 중 일부의 오류를 수정한 것이다(Le Rousseau's, 2001, 참고문헌 9 참조).
따라서 주 심볼(Principal symbol) (a)로부터 수직 느리기항
Figure PCTKR2015010951-appb-I000023
의 2차 오더 항으로 확장은 상기한 [수학식 7] 내지 [수학식 9]와 이하의 [수학식 10]에 나타낸 바와 같은 Le Rousseau's(2001)의 α[2]와
Figure PCTKR2015010951-appb-I000024
의 순환 관계식(recursive relationship)을 통해 계산할 수 있다.
[수학식 10]
Figure PCTKR2015010951-appb-I000025
여기서,
Figure PCTKR2015010951-appb-I000026
는 미소변랑의 고주파 근사(high frequency approximation) 또는 주 부분(principal part)인 3×3 행렬을 나타낸다.
Figure PCTKR2015010951-appb-I000027
와
Figure PCTKR2015010951-appb-I000028
는 이하의 [수학식 11]에 나타낸 바와 같이 대각화 행렬(diagonalizing matrix)인 M0과 이것의 역행렬인 [M0]-1의 곱으로 나타낼 수 있다.
[수학식 11]
Figure PCTKR2015010951-appb-I000029
여기서,
Figure PCTKR2015010951-appb-I000030
이고,
Figure PCTKR2015010951-appb-I000031
와
Figure PCTKR2015010951-appb-I000032
는 각각 S파와 P파의 배경 속도이고, 대각화 행렬 M0는 고유치벡터 (eigenvector)로 이루어진 행렬이며(De Hoop and De Hoop, 1994, 참고문헌 13 참조), 그 역행렬인 [M0]- 1는 얇은 막을 전파하기 직전에 배경 매질에서 P파와 S파를 분리하는 연산자로 작용한다.
즉, M0는 얇은 막 전파 후 분리된 P파와 S파를 결합하는(coupling) 역할을 하게 된다.
또한,
Figure PCTKR2015010951-appb-I000033
는 고유벡터(eigenvector) 형태로 이하의 [수학식 12]와 같이 표현된다.
[수학식 12]
Figure PCTKR2015010951-appb-I000034
[수학식 11]과 [수학식 12]로부터, 이하의 [수학식 13]이 구해진다.
[수학식 13]
Figure PCTKR2015010951-appb-I000035
따라서
Figure PCTKR2015010951-appb-I000036
는
Figure PCTKR2015010951-appb-I000037
를 구하기만 하면
Figure PCTKR2015010951-appb-I000038
를 통해 도출할 수 있다.
즉, 상기한 [수학식 11]과 [수학식 12]를 [수학식 10]에 대입하고 [수학식 8]의
Figure PCTKR2015010951-appb-I000039
를 이용하면, 이하의 [수학식 14]와 같이 하여
Figure PCTKR2015010951-appb-I000040
의 1차 확장식인
Figure PCTKR2015010951-appb-I000041
를 구할 수 있다.
마찬가지로,
Figure PCTKR2015010951-appb-I000042
의 2차 확장식
Figure PCTKR2015010951-appb-I000043
또한 상기한 [수학식 11]과 [수학식 12]를 [수학식 10]에 대입하고 [수학식 9]의
Figure PCTKR2015010951-appb-I000044
를 이용하여, 이하의 [수학식 15]와 같이 계산할 수 있다.
[수학식 14]
Figure PCTKR2015010951-appb-I000045
여기서, 나머지 성분들은 0(zero)이다.
[수학식 15]
Figure PCTKR2015010951-appb-I000046
여기서, 나머지 성분들은 0(zero)이다.
Figure PCTKR2015010951-appb-I000047
행렬의 대각 성분들은 P파와 S파의 전파에 관여하고, 나머지 성분들(off-diagonal entries)은 모드 결합(mode coupling)에 관여하게 되며, Fourier 영역에서 공간 변수 (ε)를 분리하지 않으면 모든 계산 노드에서 Fourier 변환을 실시하여 공간과 파수 영역을 오가야 하는 어려움이 있다.
그러나 EGS 전파자의 경우, 공간 변수를 효율적으로 Fourier 영역에서 분리하여 계산하기 때문에 노드의 수만큼 Fourier 변환을 하는 것이 아니라 확장된 항(term)의 수만큼만 변환을 실시하게 되어 파의 전파에 관한 계산을 빠르게 수행할 수 있다.
이와 같이 공간 변수를 주파수-파수 영역으로부터 효율적으로 분리하기 위해,
Figure PCTKR2015010951-appb-I000048
는, 이하의 [수학식 16]에 나타낸 바와 같이, 위상 항(αv)과 공간변수 항(xμ)의 곱 형태로 전개할 수 있다(Le Rousseau, 2001, 참고문헌 9 참조).
[수학식 16]
Figure PCTKR2015010951-appb-I000049
여기서, λ는 [수학식 14] 및 [수학식 15]에 있는 항(term)의 수를 나타낸다.
따라서 [수학식 11] 내지 [수학식 13] 및 [수학식 16]을 상기한 [수학식 1]과 [수학식 2]에 대입하면, 이하의 [수학식 17]과 같은 최종적인 EGS 전파자를 도출할 수 있다.
[수학식 17]
Figure PCTKR2015010951-appb-I000050
여기서,
Figure PCTKR2015010951-appb-I000051
는 스크린(screen) 항으로 주파수-공간(f-x) 영역에서 계산되고,
Figure PCTKR2015010951-appb-I000052
는 위상 이동(phase-shift) 항으로 주파수-파수(f-k) 영역에서 계산이 수행된다.
또한, N은 정규화 연산자(normalizing operator)로 무한급수 전개시 절단오차로 인해 파수 또는 전파각에 따라 에너지가 부정확하게 증폭 또는 감쇠되는 현상을 방지하기 위해 사용된다(Le Rousseau and De Hoop, 2003, 참고문헌 8 참조).
즉, 도 1을 참조하면, 도 1은 EGS 전파자의 개념을 나타내는 도면으로, 상기한 [수학식 17]의 계산을 개념적으로 나타내는 도면이다.
즉, 도 1에 나타낸 바와 같이, P파와 S파로 분리된 파동장은 각각 두 개의 배경 매질(P파와 S파의 배경 속도가 각각
Figure PCTKR2015010951-appb-I000053
인 매질과
Figure PCTKR2015010951-appb-I000054
인 매질의 얇은 막)로 진행한다.
S파가 진행하는 두 번째
Figure PCTKR2015010951-appb-I000055
매질의 경우, P파의 배경속도로
Figure PCTKR2015010951-appb-I000056
가 아닌
Figure PCTKR2015010951-appb-I000057
를 사용하는 것은 S파의 속도가 P파의 배경 속도보다 낮기 때문에 브랜치 포인트(branch point)가 발생할 가능성이 있기 때문이다.
즉, 도 1에서 있어서, P파는
Figure PCTKR2015010951-appb-I000058
의 배경속도를 가진 매질을 진행하게 되는데 이때 발생하는 S파(그림에서
Figure PCTKR2015010951-appb-I000059
)는 모드 전환에만 기여를 하고 실제 다음 얇은 막에서는 사용되지 않는다.
마찬가지로, S파는
Figure PCTKR2015010951-appb-I000060
의 배경속도를 가진 매질을 전파하는데, 이때 발생한 P파(
Figure PCTKR2015010951-appb-I000061
)는 모드 전환에만 사용이 되고 다음 막에서는 사용되지 않는다.
즉,
Figure PCTKR2015010951-appb-I000062
의 배경속도를 가진 매질을 전파한 P파와
Figure PCTKR2015010951-appb-I000063
의 배경속도를 가진 매질을 전파한 S파가 실제 다음 막으로 진행하는 P파와 S파가 된다.
또한, 상기한 바와 같이, 얇은 막을 통과한 P파와 S파 파동장은 M0에 의해 두 모드가 합쳐져 원래의 파동장(Vx, Vz)이 되고, 공간 영역(F-X 영역)에서
Figure PCTKR2015010951-appb-I000064
에 의해 스크린 항(screen term)을 계산하게 된다.
이는 다시 [M0]-1에 의해 P파와 S파로 분리되고, 다음 막에서의 입력이 e도되며, 이때 발생한 P파와 S파(그림에서 P 와 S)는 본 발명에 따른 구조보정 알고리즘에서 별도의 배열에 각 심도별 분리된 파동장으로 저장되고, 이렇게 모든 심도별로 저장된 P파와 S파는 구조보정에 있어서 영상화 조건의 입력으로 사용된다.
따라서 이러한 과정은, 구조보정 전에 입력 파동장을 별도로 P파와 S파로 분리해야하는 다른 구조보정 기법과 달리, 별도의 모드 분리 없이 바로 다성분 탄성파 자료를 구조보정의 입력으로 사용할 수 있는 핵심 단계가 된다.
계속해서, 도 2를 참조하여, 상기한 바와 같이 하여 구현된 EGS 전파자의 임펄스 응답을 검증한 내용에 대하여 설명한다.
즉, 도 2를 참조하면, 도 2는 균질한(zero-pertubation) 매질에서 본 발명의 실시예에 따른 EGS 전파자로 전파한 파동장을 나타내는 도면이다.
즉, 본 발명자들은, 상기한 [수학식 17] 및 도 1에 나타낸 바와 같은 알고리즘으로 구현한 본 발명의 실시예에 따른 EGS 전파자의 임펄스 응답(impulse response)을 검증하기 위해, P파 속도가 2,100 m/s이고 S파 속도가 1,050 m/s이며, 밀도가 2.2 g/cm3인 미소변량(perturbation)이 없는 등방성 균질 모델 상의 지표에서 수직 점 음원(point source)을 발생시켜 0.2초 후의 임펄스 반응을 나타내었다.
더 상세하게는, 도 2a는 수평 입자 속도장(horizontal particle velocity field)(Vx)의 기록이고, 도 2b는 수직 입자 속도장(vertical particle velocity field)(Vz)의 기록을 각각 나타내고 있다.
도 2a에 있어서, 수직 임펄스 음원(impulse source)일 때 수평 입자 속도 필드에서의 P파와 S파는 음원(source) 위치에서 그 극성이 바뀌게 되고, 반대로, ㄷ도 2b와 같이, 수직 입자 속도장에서는 극성이 바뀌지 않는 것을 확인할 수 있다.
또한, 탄성(elastic) 매질에서는 수직 또는 수평 점 음원(point source)이 발생하면 P파와 함께 S파도 발생이 되며, 그들의 극성은 입사각과 음원(source)의 발생 방향에 따라 달라지는데, 도 2를 참조하면, EGS 전파자로 계산한 파동장의 극성이 정확하게 묘사됨을 확인할 수 있다.
아울러, 본 발명자들은, 미소변량이 존재하는 매질에서 상기한 바와 같이 하여 구현된 본 발명의 실시예에 따른 EGS 전파자로 계산된 확장 오더(order)별 전파자를 나타내었다.
즉, 도 3을 참조하면, 도 3은 미소변량이 존재하는 매질에서 본 발명의 실시예에 따른 EGS 전파자로 전파한 파동장을 나타내는 도면이다.
더 상세하게는, 도 3a는 일정한 미소 변량이 주어진 모델 상에서 수직 점 음원이 주어졌을 때 P파의 파동장이고, 도 3b는 일정한 미소 변량이 주어진 모델 상에서 수직 점 음원이 주어졌을 때 S파의 파동장을 각각 나타내고 있다.
여기서, 각각의 P파와 S파 파동장은 상기한 알고리즘의 [M0]-1에 의해 P파와 S파로 분리된 파동장을 나타내며, 이 모델에서의 P파 속도와 S파 속도는 각각 2100 m/s와 700 m/s 이고, 밀도는 2.2 g/cm3이고, 속도와 밀도의 미소변량(perturbation)이 주어졌을 때 파형의 근사(approximation) 정도를 확인하기 위해 매질의 P파와 S파의 배경 속도는 실제 속도의 2/3로 부여하였다.
또한, 도 3에 있어서, 탄성 스플릿 스텝(Elastic Split-step)은 0차 확장을 나타내고, EGSP1은 1차 오더 확장, EGSP2는 본 발명에서 제안한 2차 오더 확장된 전파자를 각각 나타내고 있다.
즉, 도 3에 나타낸 바와 같이, 확장 오더의 차수가 클수록 보다 실제 전파에 가깝게 파가 전파함을 확인할 수 있으며, 2차 오더 전파자가 실제 파형(Exact)에 제일 가깝게 전파함을 확인할 수 있다.
다음으로, 도 4를 참조하면, 도 4는 단순 2층 수평모델에서 본 발명의 실시예에 따른 EGS 전파자로 전파한 전파시킨 각 성분(component)별 파동장과 모드 분리 연산자(operator)로 P파와 S파를 분리한 파동장을 각각 나타내는 도면이다.
즉, 본 발명자들은, EGS 전파자에 포함된 모드 분리 연산자(operator)인 [M0]-1에 의해 파동장이 P파와 S파로 정확이 분리되는지를 단순한 수평 2층 모델 상에서 확인하였다.
여기서, 사용한 모델의 상부층의 P파 속도, S파 속도 및 밀도는 각각 2500 m/s, 1200 m/s 및 2.2 g/cm3 이고, 하부층의 P파속도, S파 속도 및 밀도는 각각 3500 m/s, 1500 m/s 및 2.4 g/cm3 이며, 모델의 크기는 3×2 km, 상부층의 심도는 0.8 km 이다.
아울러, 도 4에 있어서, 지표면 1.5 km 지점에서 수평 포인트 소스(horizontal point source)를 발생시켰을 때 수직 입자 속도 필드(vertical particle velocity field)와 수평 입자 속도 필드(horizontal particle velocity field)에서의 파동장은 도 4a에 나타내었고, P파와 S파로 분리된 파동장은 도 4b에 나타내었으며, 좌측 패널의 그림은 음원이 발생한 지 0.5초 후의 파동장의 스냅샷(snap shot)이고, 오른쪽 패널의 그림은 음원이 발생한지 0.8초 후의 파동장을 나타낸다.
도 4에 나타낸 바와 같이, 수직 입자 속도 필드(Z-component)와 수평 입자 속도 필드(X-component) 상에서는 P파와 S파가 공존하며 전파해나가는 것을 확인 할 수 있고, EGS 전파자 내의 [M0]-1에 의해 모드가 분리된 단면상에서는 P파는 오직 P파 파동장에, S파는 오직 S파 파동장에만 존재함을 확인할 수 있다.
또한, 경계면에서 발생한 모드전환파, 즉, SP와 PS또한 각각 P파와 S파 파동장에만 존재하며, 즉, 본 발명의 실시예에 따른 EGS 전파자로 전파한 파동장은 각 전파 단계에서 P파와 S파가 정확히 분리됨으로써 이를 별도의 파동장 배열에 저장하여 구조보정에 사용하는 것에 의해 최종적으로 P파 구조보정 단면과 S파 구조보정 단면을 얻을 수 있다.
계속해서, 상기한 바와 같은 본 발명의 실시예에 따른 EGS 전파자를 이용하여 S파 극성 반전 보정 및 영상화 조건을 부여한 구조보정 알고리즘을 구현하는 방법에 대하여 설명한다.
일반적으로, 음원 영역에서 전파자에 의해 발생된 파동장과 수진기 영역(즉, 구조보정을 위해 사용되는 입력 다성분(multicomponent) 파동장)에서 전파자에 의해 발생된 파동장을 상호 상관(cross-correlation) 또는 컨벌루션(convolution) 시켜줌으로써 영상화가 이루어지며, 이는, 곧 최종 구조보정 단면을 얻게 됨을 의미한다.
이때, 어떠한 방식으로 상호상관 또는 컨벌루션을 시키느냐에 대한 문제가 영상화 조건(imaging condition)이며, 일반적으로, 취득한 다성분 탄성파 자료에서 P파와 S파가 정확히 분리된다면, 이하의 [수학식 18]에 나타낸 바와 같은 영상화 조건을 이용하는 것에 의해 음향 파동방정식을 이용한 구조보정 방법도 좋은 결과를 도출할 수 있다.
[수학식 18]
Figure PCTKR2015010951-appb-I000065
여기서,
Figure PCTKR2015010951-appb-I000066
는 푸리에 영역에서의 스칼라(scalar) 음원 파동장(source wavefield)을 나타내고,
Figure PCTKR2015010951-appb-I000067
은 푸리에 영역에서의 스칼라 수진기 파동장(receiver wavefield)을 나타내며, 기호 '*'는 켤레복소수(complex conjugate)를 나타낸다.
또한, 다성분 자료를 상기한 [수학식 18]에 적용한다면, 음원 파동장에서는 P파를 발생시키고 수진기 파동장에는 다성분 입력자료로부터 완벽히 분리된 S파 파동장을 입력하는 경우 최종적으로 PS 영상을 얻을 수 있다.
그러나 실제 현장에서 취득한 다성분 자료로부터 P파와 S파를 정확히 분리하기는 쉽지 않으며, 부정확하게 분리된 P파와 S파를 사용하여 상기한 [수학식 18]의 영상화 조건을 이용한다면 이미지의 왜곡을 발생시키게 된다.
반면, 본 발명의 실시예에 따른 EGS 구조보정 기법은 EGS 전파자 내의 모드 분리 연산자([M0]-1)에 의해 P파와 S파가 전파단계에서 분리되어 별도로 저장되므로, 이와 같이 분리된 P파와 S파 파동장을 상기한 [수학식 18]과 같은 영상화 조건에 대입하면 매우 양호한 결과를 도출할 수 있다.
더욱이, 모드 분리와 영상화가 동일한 주파수-공간 영역(frequency-spatial domain)에서 수행되므로, 추가적인 Fourier 변환이 필요 없어 수치적 FFT를 수행함으로 인해 발생할 수 있는 알리아싱(aliasing)이나 잡음으로부터 자유롭고 고품질의 영상을 얻을 수 있다.
이를 위해, 본 발명의 실시예에 따른 EGS 구조보정 알고리즘은, 상기한 [수학식 18]을 그대로 사용하지 않고, 보다 안정적인 영상화를 위해 Schleicher et al. 2008(참고문헌 15 참조)이 제안한 안정화된 나눔법(stabilized division method)을 적용하여, 이하의 [수학식 19]에 나타낸 바와 같은 영상화 조건을 사용하였다.
[수학식 19]
Figure PCTKR2015010951-appb-I000068
여기서, 분모항에 나타나 있는 ε는 수진기 파동장의 절대값의 제곱에 일정 스케일러(scaling fraction)(0 < λ < 1)를 적용한 값으로
Figure PCTKR2015010951-appb-I000069
로 표현할 수 있으며, i 와 j 는 음원과 수진기의 벡터 파동장(vector wavefield)(즉, 다성분 탄성파 자료에서 수평성분(x-component)과 수직성분(z-component)을 의미함)을 나타낸다.
즉, 도 5를 참조하면, 도 5는 도 1에 나타낸 EGS 전파자를 이용한 EGS 구조보정 알고리즘의 개념을 나타내는 도면이다.
더 상세하게는, 도 5에 나타낸 EGS 구조보정 알고리즘의 개념도는 도 1을 참조하여 설명한 얇은 막 전파를 전체 모델로 확장한 경우라고 할 수 있으며, 도 5에 있어서, 하향전파, 즉, "Forward Propagation"은 음원 영역에서의 파동장 전파를, 상향전파, 즉, "Backward propagation"은 수진기 영역에서의 파동장 전파를 나타낸다.
도 5에 나타낸 바와 같이, 각각의 모드 분리와 모드 결합이 각각의 얇은 막 전파 전후에 이루어지고 있고, 모드가 분리된 P파와 S파는 별도의 배열에 저장하여 최종적으로 영상화 조건에 사용하게 된다.
또한, 도 5 가운데의 기호 "*"는 영상화 조건의 상호상관을 나타내며, 이때, 영상화조건으로는 상기한 [수학식 19]가 사용되었다.
아울러, 도 6을 참조하면, 도 6은 도 5에 나타낸 바와 같은 개념도를 바탕으로 구현된 본 발명의 실시예에 따른 EGS 구조보정 알고리즘의 전체적인 구성을 개략적으로 나타내는 플로차트이다.
도 6에 나타낸 바와 같이, 본 발명의 실시예에 따른 구조보정 알고리즘은, 크게 나누어, 분석하고자 하는 탄성파 다성분 자료를 입력받아 음원(source)과 수진기(Receivers)에 대한 모델을 수립하고, 푸리에 변환을 통하여 계산해야 할 주파수 대역(frequency band)을 결정하는 단계와, 상기한 EGS 전파자를 이용하여 각각의 주파수 대역에 대하여 음원 영역에서의 파동장 전파(forward propagation from source) 및 수진기 영역에서의 파동장 전파(backward propagation from receiver)를 각각 산출하는 단계와, 상기 단계에서 각각 계산된 결과를 상호상관(Cross Correlation)을 통해 통합하고, 상기한 영상화 조건(Imaging condition)에 의해 구조보정된 영상 데이터를 출력하는 일련의 단계를 포함하여 구성될 수 있다.
여기서, 상기한 음원 영역에서의 파동장 전파를 산출하는 단계는, 상기한 본 발명의 실시예에 따른 EGS 전파자 내의 모드 분리 연산자(operator)([M0]-1)에 의해 음원 파동장(source wavefield)을 분리하는 단계와, f-x 영역에서 스크린(screen) 및 모드 결합(mode coupling)을 계산하는 단계와, 계산된 결과에 공간 변수(spatial variable, x)에 대하여 푸리에 변환을 수행한 후, f-k 영역에서 외삽된(extrapolated) 파동장을 계산하는 단계와, 구조보정(migration)을 위해 각각의 모드를 저장하고, 상기한 EGS 전파자 내의 모드결합 연산자(M0)에 의해 음원 파동장을 재구성(recomposition)하는 단계 및 재구성된 파동장에 역푸리에 변환을 수행하는 단계를 포함하여 구성될 수 있다.
또한, 상기한 수진기 영역에서의 파동장 전파를 산출하는 단계는, 마찬가지로, 상기한 모드 분리 연산자(operator)([M0]-1)에 의해 수진기 파동장(receiver wavefield)을 분리하는 단계와, f-x 영역에서 스크린(screen) 및 모드 결합(mode coupling)을 계산하는 단계와, 계산된 결과에 공간 변수(spatial variable, x)에 대하여 푸리에 변환을 수행한 후, f-k 영역에서 외삽된(extrapolated) 파동장을 계산하는 단계와, 구조보정(migration)을 위해 각각의 모드를 저장하고, 상기한 EGS 전파자 내의 모드결합 연산자(M0)에 의해 수진기 파동장을 재구성(recomposition)하는 단계 및 재구성된 파동장에 역푸리에 변환을 수행하는 단계를 포함하여 구성될 수 있다.
아울러, 상기한 EGS 전파자는, 상기한 [수학식 1] 내지 [수학식 17]을 참조하여 설명한 바와 같이 하여 구해진 EGS 전파자를 이용할 수 있다.
더욱이, 상기한 영상화 조건(Imaging condition)은, 상기한 [수학식 18] 내지 [수학식 19]를 참조하여 설명한 바와 같은 영상화 조건을 이용할 수 있다.
또한, 도 6에 있어서, 각각의 코드는 병렬처리(parallel processing)가 가능하도록 MPI(Message Passing Interface) 코드로 최적화하였으며, 이때, MPI를 위한 플로차트는 도 7에 나타낸 바와 같이 하여 구성될 수 있다.
즉, 도 7을 참조하면, 도 7은 도 6에 나타낸 본 발명의 실시예에 따른 EGS 구조보정 알고리즘의 병렬처리를 위한 MPI를 구현하는 과정의 전체적인 구성을 개략적으로 나타내는 플로차트이다.
도 7에 나타낸 바와 같이, 복수의 프로세서를 이용하여, 각각의 주파수 대역별로 처리를 할당하여 병렬처리가 이루어지도록 구성함으로써, 전체적인 처리시간을 단축할 수 있다.
여기서, S파는 매질을 전파하면서 매질 경계면의 반사점에서 그 극성이 변하게 되며, 종래의 탄성(elastic) 구조보정 기법들은 대부분 이러한 S파 극성 역전 현상을 공간 영역(spatial domain)에서 분리해야 하나, 반사면이 수평이 아닌 복잡한 실제 지하구조의 경우, 반사 지점이 정확히 어디인지 규명하기가 쉽지 않다는 문제점이 있었다.
반면, 본 발명의 실시예에 따른 EGS 기법은, 전파자 계산의 특성상 주파수-파수 영역에서 그 극성 역전 현상에 대한 보정을 실시할 수 있다.
더 상세하게는, Sava and Fomel. 2003(참고문헌 14 참조)에 의하면, 파수 영역에서 반사 지점에서의 반사각은 이하의 [수학식 20]에 의해 구할 수 있다.
[수학식 20]
Figure PCTKR2015010951-appb-I000070
여기서, γ는 반사지점에서의 반사각을 나타내며 kh와 kz는 각각 거리 방향의 파수와 심도 방향의 파수를 나타낸다.
상기한 [수학식 20]에 의해, 반사지점에서의 반사각의 부호(sign)는 -1×수평 파수의 부호(sign of horizontal wavenumber)로 표현할 수 있고, 따라서 S파의 극성은 주파수-파수 영역에서 kx = 0인 지점을 기준으로 한쪽의 부호를 반대 부호로 바꾸어 주기만 하면 된다.
즉, 도 8을 참조하면, 도 8은 상기한 [수학식 20]에 나타낸 영상화 조건을 이용하여 극성 역전의 보정 방법을 적용한 경우와 적용하지 않은 경우의 구조보정 결과를 각각 비교하여 나타내는 도면이다.
도 8에 있어서, 도 8a는 단순 2층 모델에서 한 개의 음원으로 발생시킨 합성 탄성파 자료(synthetic seismic data)를 이용하여 본 발명의 실시예에 따른 EGS 구조보정 기법을 통해 구조보정을 실시한 결과이고, 도 8b는 PP 단면, 도 8c는 극성 역전 현상을 보정하지 않은 PS단면, 도 8d는 상기한 [수학식 20]을 이용하여 극성 역전현상을 보정한 PS 단면을 각각 나타내고 있다.
따라서 도 8에 나타낸 바와 같이, 도 8c에는 반사점 지점에서 PS의 극성이 바뀜을 확인할 수 있고, 그 보정 단면인 도 8d에는 극성 역전현상이 없음을 확인할 수 있다.
계속해서, 상기한 바와 같이 하여 구현되는 본 발명의 실시예에 따른 EGS 구조보정 방법을 실제 모델에 적용하여 검증을 실시한 결과에 대하여 설명한다.
즉, 본 발명자들은, 상기한 EGS 구조보정 기법을 보다 현실적이고 복잡한 모델에 적용하여 실제 모델과 구조보정을 통해 도출한 영상이 잘 맞는지 확인하였으며, 이를 위해, 탄성파 탐사 분야에서 벤치마크(Benchmark) 모델로 널리 알려진 SEG/EAGE Salt Model(Aminzadeh et al., 1997, 참고문헌 16 참조)에서 수치 모델링을 통해 인공적으로 다성분 탄성파 자료(synthetic data)를 생성하고, 이를 이용하여 본 발명의 실시예에 따른 EGS 구조보정 기법을 적용하여 구조보정 단면을 얻었다.
여기서, SEG/EAGE Salt 모델은 본래 3차원 모델이나, 본 실시예에서는 x = 7680m인 부분의 2차원 단면 모델을 추출하고 2차원 모델링을 수행하여 자료를 생성하였다.
즉, 도 9를 참조하면, 도 9는 상기한 바와 같이 하여 생성된 2차원 단면 모델의 P파 속도 구조를 나타내는 도면이다.
도 9에 있어서, S파 속도 구조는 P파 속도에
Figure PCTKR2015010951-appb-I000071
을 일괄적으로 곱하여 생성하였으며, 밀도 모델은 Gardner 식을 이용하여 도출하였다.
또한, 이 모델의 지표상의 760 ~ 12760m 영역에 40m 간격으로 음원을 설치하고, 수진기는 60 ~ 13460m 영역에 20m 간격으로 설치하여 합성 탄성파 자료를 취득하였으며, 이때, 음원의 주 주파수는 5Hz이고 최대 주파수는 12.5Hz이다.
아울러, 구조보정을 위한 그리드(grid) 간격은 S파 최소속도를 고려하여 dx = 20m, dz = 10m로 설정하고, 0 ~ 12Hz까지의 주파수 영역으로 구조보정을 실시하여 최종 PP 및 PS 단면을 도출하였다.
즉, 도 10을 참조하면, 도 10은 본 발명의 실시예에 따른 EGS 구조보정 기법을 적용하여 도출한 최종 PP 및 PS 단면의 영상을 나타내는 도면이다.
도 10에 있어서, 먼저, 도 10a는 PP 구조보정 결과이며, 이는 도 9의 속도 영상과 비교했을 때 암염돔(Salt dome)의 영상이 매우 정확하게 표현되었고, 주변 반사 경계면도 정확한 위치에 영상화되어 있음을 확인할 수 있다.
또한, 도 10b는 S파 극성 역전 현상을 보정하지 않고 구조보정을 실시한 PS 결과이고 도 10c는 본 발명의 실시예에 따른 EGS 구조보정 알고리즘에 포함된 극성 역전 보정 모듈을 사용한 PS 결과이다.
도 10b에 나타낸 바와 같이, 극성 역전 현상을 보정하지 않은 경우 이미지의 품질이 떨어지고, 특히, 암염(Salt) 상부 이벤트의 경우 매끈한 곡선이 아닌 점 형태로 나타나 있으며, 또한, 암염(Salt) 하부 가로로 긴 수평 이벤트의 경우 극성 역전 보정을 하지 않은 결과에는 잘 나타나지 않는다.
즉, S파 극성 역전 현상을 보정하지 않은 경우는 도 9에 나타낸 원래의 속도 모델과 비교했을 때 이미지의 연속성이 떨어지고 암염(Salt) 하부의 수평층은 아예 보이지도 않는데, 이는 각 음원 모음(shot gather) 별로 구조보정하고 추후 적분(또는 합)하는 과정에서 극성이 역전하는 부근에서 서로 상쇄되어 이미지의 품질을 저하시킨 결과이다.
반면, 도 10c에 나타낸 바와 같이, 본 발명의 실시예에 따른 EGS 구조보정 알고리즘에 포함된 극성 역전 보정 모듈을 통해 S파 극성 역전 현상을 보정하여 구조보정을 실시한 결과는, 도 9의 원래 속도 모델이나 도 10a의 PP결과와 비교했을 때 정확한 암염(Salt body)을 영상화하였고, 주변의 반사 이벤트들도 비교적 정확한 위치에 영상화되었음을 확인할 수 있다.
상기한 바와 같은 내용으로부터, 본 발명의 실시예에 따른 EGS 전파자 및 이를 이용한 EGS 구조보정 방법은 종래기술의 방법들에 비해 더욱 정확하고 사실적인 구조보정 및 영상화가 가능하여 구조보정 영상의 품질을 향상시킬 수 있는 것임을 알 수 있다.
따라서 상기한 바와 같이 하여 본 발명에 따른 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법을 구현할 수 있다.
또한, 상기한 바와 같이 하여 본 발명에 따른 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법을 구현하는 것에 의해, 본 발명에 따르면, 수직 느리기 항(vertical slowness operator)을 2차항까지 확장하여 복잡한 모델이나 수평적으로 속도 변화가 심한 모델에서도 보다 정확한 파의 전파를 구현하여 EGS 전파자의 성능을 더욱 향상시킬 수 있도록 구성되는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법이 제공됨으로써, 수직 느리기 항이 1차항 까지만 근사되는 한계가 있었던 종래기술의 EGS 전파자의 문제점을 해결할 수 있다.
또한, 본 발명에 따르면, EGS 전파자 내부에 구현된 P-S분리 모듈을 포함하여 입력 다성분 파동장을 P파와 S파로 분리할 필요 없이 음원 모음(shot gather)을 바로 구조보정 입력으로 사용하여 P파 및 S파 영상 단면을 생성할 수 있도록 구성되는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법이 제공됨으로써, 구조보정 전에 파동장을 P파와 S파로 분리하여 입력으로 사용함으로 인해 계산 및 구조가 복잡해지는 단점이 있었던 종래기술의 구조보정 기법들의 문제점을 해결할 수 있다.
아울러, 본 발명에 따르면, S파 영상화 전에 주파수-파수 영역에서 극성 전환에 대한 보정을 수행하는 모듈을 추가하여 S파 구조보정 영상의 품질을 향상시킬 수 있도록 구성되는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법을 제공됨으로써, 탄성파 S파 음원 모음(shot gather)의 경우 반사지점에서 극성이 역전되는 현상이 있음으로 인해 이를 그대로 구조보정에 사용하게 되면 이벤트의 연속성이 떨어지게 되어 최종 단면의 품질이 저하되는 문제가 있었던 종래기술의 구조보정 기법들의 문제점을 해결할 수 있다.
이상, 상기한 바와 같은 본 발명의 실시예를 통하여 본 발명에 따른 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법의 상세한 내용에 대하여 설명하였으나, 본 발명은 상기한 실시예에 기재된 내용으로만 한정되는 것은 아니며, 따라서 본 발명은, 본 발명이 속하는 기술분야에서 통상의 지식을 가진 자에 의해 설계상의 필요 및 기타 다양한 요인에 따라 여러 가지 수정, 변경, 결합 및 대체 등이 가능한 것임은 당연한 일이라 하겠다.

Claims (8)

  1. 지하 매질의 경계면들을 전파하면서 P파와 S파 상호간의 모드 전환을 거치는 탄성파의 거동을 표현하기 위한 EGS(Elastic Generalized-Screen) 전파자(wave propagator)의 수직 느리기항(vertical slowness term)을 2차항까지 확장하여 다성분 자료의 구조보정시 입력자료를 P파와 S파 파동장으로 분리할 필요 없이 음원 모음(shot gather)을 바로 구조보정 입력으로 사용하여 P파 및 S파 영상 단면을 생성할 수 있도록 구성되는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법에 있어서,
    분석하고자 하는 탄성파 다성분 자료를 입력받아 음원(source)과 수진기(receivers)에 대한 모델을 수립하고, 푸리에 변환(Fourier transform)을 통하여 계산해야 할 주파수 대역(frequency band)을 결정하는 전처리 단계;
    상기 EGS 전파자를 이용하여 각각의 주파수 대역에 대하여 음원 영역에서의 파동장 전파(forward propagation from source)를 산출하는 음원영역 처리단계;
    상기 EGS 전파자를 이용하여 각각의 주파수 대역에 대하여 수진기 영역에서의 파동장 전파(backward propagation from receiver)를 산출하는 수진기영역 처리단계;
    상기 음원영역 처리단계 및 상기 수진기영역 처리단계에서 각각 계산된 결과를 상호상관(Cross Correlation)을 통해 통합하고 영상화 조건(Imaging condition)에 의해 구조보정을 행하는 구조보정단계; 및
    상기 구조보정단계에서 구조보정된 영상 데이터를 출력하는 출력단계를 포함하여 구성되는 것을 특징으로 하는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법.
  2. 제 1항에 있어서,
    상기 EGS 전파자는,
    이하의 수학식에 나타낸 EGS 전파자를 이용하도록 구성되는 것을 특징으로 하는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법.
    Figure PCTKR2015010951-appb-I000072
    (여기서, xμ(μ = 1, 2)과 x3는 각각 수평과 수직 좌표를 나타내고, s = -iω(ω는 각주파수) 이며, αv = kv/iω(kv는 수평성분 파수, v = 1, 2)이고,
    Figure PCTKR2015010951-appb-I000073
    와
    Figure PCTKR2015010951-appb-I000074
    는 각각 Fourier 변환과 그의 역변환을 나타내며, M0는 고유치벡터 (eigenvector)로 이루어진 대각화 행렬로서 분리된 P파와 S파를 결합하는(coupling) 연산자로 작용하고, [M0]-1은 M0의 역행렬로서 P파와 S파를 분리하는 연산자로 작용하며, λ는 항(term)의 수를 나타냄.)
  3. 제 2항에 있어서,
    상기 음원영역 처리단계는,
    상기 EGS 전파자 내의 모드 분리 연산자(operator)([M0]-1)에 의해 음원 파동장(source wavefield)을 분리하는 단계;
    주파수-공간(f-x) 영역에서 스크린(screen) 및 모드 결합(mode coupling)을 계산하는 단계;
    계산된 결과에 공간 변수(spatial variable, x)에 대하여 푸리에 변환을 수행하는 단계;
    주파수-파수(f-k) 영역에서 외삽된(extrapolated) 파동장을 계산하는 단계;
    구조보정(migration)을 위해 각각의 모드를 저장하고, 상기 EGS 전파자 내의 모드결합 연산자(M0)에 의해 상기 음원 파동장을 재구성(recomposition)하는 단계; 및
    재구성된 상기 음원 파동장에 역푸리에 변환을 수행하는 단계를 포함하여 구성되는 것을 특징으로 하는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법.
  4. 제 3항에 있어서,
    상기 수진기영역 처리단계는,
    상기 EGS 전파자 내의 모드 분리 연산자(operator)([M0]-1)에 의해 수진기 파동장(receiver wavefield)을 분리하는 단계;
    주파수-공간(f-x) 영역에서 스크린(screen) 및 모드 결합(mode coupling)을 계산하는 단계;
    계산된 결과에 공간 변수(spatial variable, x)에 대하여 푸리에 변환을 수행하는 단계;
    주파수-파수(f-k) 영역에서 외삽된(extrapolated) 파동장을 계산하는 단계;
    구조보정(migration)을 위해 각각의 모드를 저장하고, 상기 EGS 전파자 내의 모드결합 연산자(M0)에 의해 상기 수진기 파동장을 재구성(recomposition)하는 단계; 및
    재구성된 상기 수진기 파동장에 역푸리에 변환을 수행하는 단계를 포함하여 구성되는 것을 특징으로 하는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법.
  5. 제 4항에 있어서,
    상기 음원영역 처리단계 및 상기 수진기영역 처리단계는,
    복수의 프로세서를 이용하여, 각각의 주파수 대역별로 처리를 할당하여 병렬처리(parallel processing)가 가능하도록 구성됨으로써, 전체적인 처리시간을 단축할 수 있도록 구성되는 것을 특징으로 하는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법.
  6. 제 5항에 있어서,
    상기 구조보정단계는,
    상기 영상화 조건으로서, 이하의 수학식에 나타낸 영상화 조건을 이용하도록 구성되는 것을 특징으로 하는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법.
    Figure PCTKR2015010951-appb-I000075
    (여기서, Iij는 최종 영상화 이미지,
    Figure PCTKR2015010951-appb-I000076
    는 푸리에 영역에서의 스칼라(scalar) 음원 파동장(source wavefield),
    Figure PCTKR2015010951-appb-I000077
    은 푸리에 영역에서의 스칼라 수진기 파동장(receiver wavefield)을 나타내고, ε는
    Figure PCTKR2015010951-appb-I000078
    로 표현할 수 있으며, i 와 j 는 음원과 수진기의 벡터 파동장(vector wavefield)으로, 다성분 탄성파 자료에서 수평성분(x-component)과 수직성분(z-component)을 각각 의미함.)
  7. 제 6항에 있어서,
    상기 방법은,
    이하의 수학식을 이용하여 주파수-파수 영역에서 반사 지점에서의 반사각을 구하는 것에 의해, 매질 경계면의 반사점에서 극성이 변화하게 되는 S파 극성 역전 현상을 보정하는 극성보정단계를 더 포함하여 구성되는 것을 특징으로 하는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법.
    Figure PCTKR2015010951-appb-I000079
    (여기서, γ는 반사지점에서의 반사각을 나타내며, kh와 kz는 각각 거리 방향의 파수와 심도 방향의 파수를 나타냄.)
  8. 청구항 1항 내지 청구항 7항 중 어느 한 항에 기재된 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법을 이용하여, 다성분 자료를 P파와 S파로 분리할 필요 없이 음원 모음(shot gather)을 바로 구조보정 입력으로 사용하여 P파 및 S파 영상 단면을 생성할 수 있도록 구성되는 동시에, S파 영상화 전에 파수-주파수 영역에서 극성 전환에 대한 보정을 수행하여 S파 구조보정 영상의 품질을 향상시킬 수 있도록 구성되는 것을 특징으로 하는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 시스템.
PCT/KR2015/010951 2014-10-17 2015-10-16 탄성파 다성분 자료에 대한 중합전 egs 구조보정 방법 WO2016060513A1 (ko)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US15/508,957 US20170299745A1 (en) 2014-10-17 2015-10-16 Prestack egs migration method for seismic wave multi-component data

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
KR1020140140533A KR101549388B1 (ko) 2014-10-17 2014-10-17 탄성파 다성분 자료에 대한 중합전 ｅｇｓ 구조보정 방법
KR10-2014-0140533 2014-10-17

Publications (1)

Publication Number Publication Date
WO2016060513A1 true WO2016060513A1 (ko) 2016-04-21

Family

ID=54246978

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/KR2015/010951 WO2016060513A1 (ko) 2014-10-17 2015-10-16 탄성파 다성분 자료에 대한 중합전 egs 구조보정 방법

Country Status (3)

Country Link
US (1) US20170299745A1 (ko)
KR (1) KR101549388B1 (ko)
WO (1) WO2016060513A1 (ko)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105938203A (zh) * 2016-06-24 2016-09-14 中国石油天然气股份有限公司 一种储层特性的检测方法及装置
CN106772585A (zh) * 2017-01-26 2017-05-31 中国科学院地质与地球物理研究所 一种基于弹性波解耦方程的优化拟解析方法及装置
CN106896403A (zh) * 2017-05-05 2017-06-27 中国石油集团东方地球物理勘探有限责任公司 弹性高斯束偏移成像方法和系统
CN108072896A (zh) * 2016-11-18 2018-05-25 中国石油化工股份有限公司 一种全自动地震波初至拾取方法及系统

Families Citing this family (24)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105629301B (zh) * 2016-03-29 2018-02-09 中国地质大学（北京） 薄层弹性波反射系数快速求解方法
KR101907635B1 (ko) 2016-05-30 2018-10-12 한양대학교 산학협력단 역시간 구조보정을 이용하는 지하구조 영상장치 및 방법
CN107918156B (zh) * 2017-10-30 2019-09-06 中国石油天然气集团公司 检测海底节点采集地震数据极性的方法及装置
KR101865016B1 (ko) * 2018-01-16 2018-06-05 채휘영 지표투과레이다 자료 기반의 공동 검출 방법 및 공동 검출 장치
CN110058300B (zh) * 2018-10-30 2022-06-17 南方科技大学 一次波重建方法、装置、终端设备及存储介质
CN111812708B (zh) * 2019-04-11 2022-08-30 中国石油天然气股份有限公司 地震波成像方法及装置
CN111488638B (zh) * 2020-03-13 2024-04-05 天津大学 周期分布排桩对平面sv波散射解析解的求解方法
CN112782761A (zh) * 2020-10-16 2021-05-11 中国石油大学(华东) 单程波正演模拟方法、装置、存储介质及处理器
CN112363224A (zh) * 2020-10-29 2021-02-12 中国石油天然气集团有限公司 一种横波地震勘探的应用潜力判断方法及装置
CN112462426B (zh) * 2020-11-02 2024-05-28 中国石油天然气集团有限公司 横波矢量静校正方法及装置
CN112379430B (zh) * 2020-11-13 2024-02-13 中国地质科学院 一种角度域多分量偏移成像方法
CN112462428B (zh) * 2020-11-13 2024-02-20 中国地质科学院 一种多分量地震资料偏移成像方法及系统
CN112698400B (zh) * 2020-12-04 2023-06-23 中国科学院深圳先进技术研究院 反演方法、反演装置、计算机设备和计算机可读存储介质
US20220283329A1 (en) * 2021-03-05 2022-09-08 Aramco Overseas Company B.V. Method and system for faster seismic imaging using machine learning
CN113030266A (zh) * 2021-03-25 2021-06-25 武汉理工大学 汽车三代轮毂轴承外圈超声相控阵检测装置及方法
US11940585B2 (en) 2021-04-06 2024-03-26 Saudi Arabian Oil Company System and method for estimating one-way propagation operators
CN113156514B (zh) * 2021-04-25 2022-08-23 中南大学 基于主频波数域均值滤波的地震数据去噪方法及系统
CN113406698B (zh) * 2021-05-24 2023-04-21 中国石油大学(华东) 一种基于纵横波解耦的双相介质弹性波逆时偏移成像方法
CN113933900A (zh) * 2021-10-15 2022-01-14 张焕钧 一种隧道超前探测成像方法及装置
CN116413801B (zh) * 2023-02-13 2023-09-29 中国石油大学(北京) 一种各向异性介质弹性波高精度成像方法和系统
CN115993650B (zh) * 2023-03-22 2023-06-06 中国石油大学(华东) 一种基于棱柱波的地震干涉成像方法
CN117233838B (zh) * 2023-09-20 2024-04-05 长江大学 一种二维vti介质中的弹性准纵横波场分离和逆时偏移成像方法
CN117518265A (zh) * 2023-11-01 2024-02-06 大连理工大学 层状介质场地下基于频散特性的多模态面波波场反演方法
CN117249894B (zh) * 2023-11-16 2024-04-05 自然资源部第一海洋研究所 一种水下远场声传播在海底透射厚度的诊断方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20100029214A (ko) * 2007-06-26 2010-03-16 신창수 지하 구조 영상화를 위해 라플라스 도메인에서의 파형 역산을 이용한 속도 분석 방법
WO2011100077A1 (en) * 2010-02-10 2011-08-18 Exxonmobil Upstream Research Company Methods for subsurface parameter estimation in full wavefield inversion and reverse-time migration
KR101355106B1 (ko) * 2011-06-27 2014-01-22 서울대학교산학협력단 복소주파수 그룹을 기초로 하는 지하구조 영상화 방법
KR20140021584A (ko) * 2011-03-30 2014-02-20 엑손모빌 업스트림 리서치 캄파니 스펙트럼 형상을 이용하는 전 파동장 역산의 수렴 레이트

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20100029214A (ko) * 2007-06-26 2010-03-16 신창수 지하 구조 영상화를 위해 라플라스 도메인에서의 파형 역산을 이용한 속도 분석 방법
WO2011100077A1 (en) * 2010-02-10 2011-08-18 Exxonmobil Upstream Research Company Methods for subsurface parameter estimation in full wavefield inversion and reverse-time migration
KR20140021584A (ko) * 2011-03-30 2014-02-20 엑손모빌 업스트림 리서치 캄파니 스펙트럼 형상을 이용하는 전 파동장 역산의 수렴 레이트
KR101355106B1 (ko) * 2011-06-27 2014-01-22 서울대학교산학협력단 복소주파수 그룹을 기초로 하는 지하구조 영상화 방법

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
KIM, BYEONG YEOP.: "Prestack Elastic Generalized-Screen Migration for Multicomponent Data", DOCTORAL THESIS., August 2012 (2012-08-01) *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105938203A (zh) * 2016-06-24 2016-09-14 中国石油天然气股份有限公司 一种储层特性的检测方法及装置
CN108072896A (zh) * 2016-11-18 2018-05-25 中国石油化工股份有限公司 一种全自动地震波初至拾取方法及系统
CN108072896B (zh) * 2016-11-18 2019-08-27 中国石油化工股份有限公司 一种全自动地震波初至拾取方法及系统
CN106772585A (zh) * 2017-01-26 2017-05-31 中国科学院地质与地球物理研究所 一种基于弹性波解耦方程的优化拟解析方法及装置
CN106772585B (zh) * 2017-01-26 2018-11-09 中国科学院地质与地球物理研究所 一种基于弹性波解耦方程的优化拟解析方法及装置
CN106896403A (zh) * 2017-05-05 2017-06-27 中国石油集团东方地球物理勘探有限责任公司 弹性高斯束偏移成像方法和系统
CN106896403B (zh) * 2017-05-05 2019-09-13 中国石油天然气集团有限公司 弹性高斯束偏移成像方法和系统

Also Published As

Publication number Publication date
KR101549388B1 (ko) 2015-09-02
US20170299745A1 (en) 2017-10-19

Similar Documents

Publication Publication Date Title
WO2016060513A1 (ko) 탄성파 다성분 자료에 대한 중합전 egs 구조보정 방법
Wang et al. Premigration deghosting for marine streamer data using a bootstrap approach in tau-p domain
Luo et al. A deconvolution-based objective function for wave-equation inversion
Wang et al. Full-waveform inversion with the reconstructed wavefield method
Sun et al. Phase correction in separating P-and S-waves in elastic data
US11378704B2 (en) Method for generating an image of a subsurface of an area of interest from seismic data
US10036818B2 (en) Accelerating full wavefield inversion with nonstationary point-spread functions
SG182262A1 (en) Methods for subsurface parameter estimation in full wavefield inversion and reverse-time migration
Vallée Stabilizing the empirical Green function analysis: Development of the projected Landweber method
US5530679A (en) Method for migrating seismic data
US20170276814A1 (en) Method of operating a data-processing system for the simulation of the acoustic wave propagation in the transversely isotropic media comprising an hydrocarbon reservoir
Guo et al. Target-oriented waveform redatuming and high-resolution inversion: Role of the overburden
Song et al. An efficient wavefield inversion for transversely isotropic media with a vertical axis of symmetry
Feng et al. Multiscale phase inversion for vertical transverse isotropic media
Zhang et al. Optimized Chebyshev Fourier migration: A wide-angle dual-domain method for media with strong velocity contrasts
Chavent Data space reflectivity and the migration based travel time approach to FWI
Yang et al. Image-domain waveform tomography with two-way wave-equation
Qin et al. Velocity model building from waveform tomography of band limited reflection seismic data
WO2022065981A1 (ko) 동영상 처리 장치 및 방법
Liu et al. Extended reflection waveform inversion via differential semblance optimization
Xu* et al. Pure quasi-P wave calculation in anisotropic media
Nouibat et al. Ambient‐Noise Wave‐Equation Tomography of the Alps and Ligurian‐Provence Basin
Ramos-Martínez et al. Full-waveform inversion by pseudo-analytic extrapolation
Chauris et al. Image-domain versus data-domain velocity analysis based on true-amplitude subsurface extended migration
Liu et al. Linearized extended waveform inversion via differential semblance optimization in the depth-oriented extension

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 15850834

Country of ref document: EP

Kind code of ref document: A1

WWE Wipo information: entry into national phase

Ref document number: 15508957

Country of ref document: US

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 15850834

Country of ref document: EP

Kind code of ref document: A1