KR101549388B1 - 탄성파 다성분 자료에 대한 중합전 ｅｇｓ 구조보정 방법 - Google Patents

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

Info

Publication number
KR101549388B1
KR101549388B1 KR1020140140533A KR20140140533A KR101549388B1 KR 101549388 B1 KR101549388 B1 KR 101549388B1 KR 1020140140533 A KR1020140140533 A KR 1020140140533A KR 20140140533 A KR20140140533 A KR 20140140533A KR 101549388 B1 KR101549388 B1 KR 101549388B1
Authority
KR
South Korea
Prior art keywords
wave
egs
propagator
correction
elastic
Prior art date
Application number
KR1020140140533A
Other languages
English (en)
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 KR1020140140533A priority Critical patent/KR101549388B1/ko
Application granted granted Critical
Publication of KR101549388B1 publication Critical patent/KR101549388B1/ko
Priority to PCT/KR2015/010951 priority patent/WO2016060513A1/ko
Priority to US15/508,957 priority patent/US20170299745A1/en

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

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

탄성파 다성분 자료에 대한 중합전 ＥＧＳ 구조보정 방법{Prestack elastic generalized-screen migration method for seismic multicomponent data}
본 발명은 탄성파 구조보정(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 112014098966119-pat00001

(여기서, xμ(μ = 1, 2)과 x3는 각각 수평과 수직 좌표를 나타내고, s = -iω(ω는 각주파수) 이며, αv = kv/iω(kv는 수평성분 파수, v = 1, 2)이고,
Figure 112014098966119-pat00002
와
Figure 112014098966119-pat00003
는 각각 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 112014098966119-pat00004

(여기서, Iij는 최종 영상화 이미지,
Figure 112014098966119-pat00005
는 푸리에 영역에서의 스칼라(scalar) 음원 파동장(source wavefield),
Figure 112014098966119-pat00006
은 푸리에 영역에서의 스칼라 수진기 파동장(receiver wavefield)을 나타내고, ε는
Figure 112014098966119-pat00007
로 표현할 수 있으며, i 와 j 는 음원과 수진기의 벡터 파동장(vector wavefield)으로, 다성분 탄성파 자료에서 수평성분(x-component)과 수직성분(z-component)을 각각 의미함.)
또한, 상기 방법은, 이하의 수학식을 이용하여 주파수-파수 영역에서 반사 지점에서의 반사각을 구하는 것에 의해, 매질 경계면의 반사점에서 극성이 변화하게 되는 S파 극성 역전 현상을 보정하는 극성보정단계를 더 포함하여 구성되는 것을 특징으로 한다.
Figure 112014098966119-pat00008

(여기서, γ는 반사지점에서의 반사각을 나타내며, 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 112014098966119-pat00009

여기서, xμ(μ = 1, 2)과 x3는 각각 수평과 수직 좌표를 나타내고, s = -iω(ω는 각주파수) 이며, αv = kv/iω(kv는 수평성분 파수, v = 1, 2)이다.
의사 미분 연산자로 표현된 [수학식 1]은 주파수-파수 영역에서 매질의 수직느리기항(γ)을 이용하여 위상을 이동시킴으로써 심도 방향으로 파동장을 외삽 (extrapolation) 해주는 연산자로 공간 변수(spatial variable, x)와 Fourier 변수 (αv)를 하나의 항으로 표현하게 된다.
이때,
Figure 112014098966119-pat00010
와
Figure 112014098966119-pat00011
는 각각 Fourier 변환과 그의 역 변환을 나타내는 커널이다.
또한, 수평적으로 속도변화가 있는 매질에서의 수직 느리기 항은 이하의 [수학식 2]와 같이 배경 매질에서의 위상 이동만을 고려한 수직 느리기항
Figure 112014098966119-pat00012
과 수평적 속도 변화를 고려한 미소 변량(perturbation)항
Figure 112014098966119-pat00013
으로 표현된다.
[수학식 2]
Figure 112014098966119-pat00014

여기서, 상기한 [수학식 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 112014098966119-pat00015

여기서, A의 심볼인 ar을 첫째 항까지만 근사하면(high-frequency approximation)(De Hoop, 1996, 참고문헌 12 참조), 이하의 [수학식 4]에 나타낸 바와 같다.
[수학식 4]
Figure 112014098966119-pat00016

여기서, 아래 첨자는 심볼 오더(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 112014098966119-pat00017

여기서, α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 112014098966119-pat00018
의 미소변량(perturbation)은 이하의 [수학식 6]과 같이 표현할 수 있다.
[수학식 6]
Figure 112014098966119-pat00019

또한, 주 심볼 행렬(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 112014098966119-pat00020

[수학식 8]
Figure 112014098966119-pat00021

[수학식 9]
Figure 112014098966119-pat00022

여기서, 상기한 [수학식 8]은 [수학식 5]의 테일러 전개 후 1차 입실론항(first epsilon oredr term)에 의해 재정리된(rearranged) 것이고, [수학식 9]는 [수학식 5]의 테일러 전개 후 2차 입실론항에 의해 재정리된(rearranged) 것이며, 이때, 특정 항들은 루소의 유도식(Le Rousseau's derivation) 중 일부의 오류를 수정한 것이다(Le Rousseau's, 2001, 참고문헌 9 참조).
따라서 주 심볼(Principal symbol) (a)로부터 수직 느리기항
Figure 112014098966119-pat00023
의 2차 오더 항으로 확장은 상기한 [수학식 7] 내지 [수학식 9]와 이하의 [수학식 10]에 나타낸 바와 같은 Le Rousseau's(2001)의 α[2]와
Figure 112014098966119-pat00024
의 순환 관계식(recursive relationship)을 통해 계산할 수 있다.
[수학식 10]
Figure 112014098966119-pat00025

여기서,
Figure 112014098966119-pat00026
는 미소변랑의 고주파 근사(high frequency approximation) 또는 주 부분(principal part)인 3×3 행렬을 나타낸다.
Figure 112014098966119-pat00027
와
Figure 112014098966119-pat00028
는 이하의 [수학식 11]에 나타낸 바와 같이 대각화 행렬(diagonalizing matrix)인 M0과 이것의 역행렬인 [M0]-1의 곱으로 나타낼 수 있다.
[수학식 11]
Figure 112014098966119-pat00029

여기서,
Figure 112014098966119-pat00030
이고,
Figure 112014098966119-pat00031
와
Figure 112014098966119-pat00032
는 각각 S파와 P파의 배경 속도이고, 대각화 행렬 M0는 고유치벡터 (eigenvector)로 이루어진 행렬이며(De Hoop and De Hoop, 1994, 참고문헌 13 참조), 그 역행렬인 [M0]- 1는 얇은 막을 전파하기 직전에 배경 매질에서 P파와 S파를 분리하는 연산자로 작용한다.
즉, M0는 얇은 막 전파 후 분리된 P파와 S파를 결합하는(coupling) 역할을 하게 된다.
또한,
Figure 112014098966119-pat00033
는 고유벡터(eigenvector) 형태로 이하의 [수학식 12]와 같이 표현된다.
[수학식 12]
Figure 112014098966119-pat00034

[수학식 11]과 [수학식 12]로부터, 는 이하의 [수학식 13]과 같이 하여 구해진다.
[수학식 13]
Figure 112014098966119-pat00035

따라서
Figure 112014098966119-pat00036
는
Figure 112014098966119-pat00037
를 구하기만 하면
Figure 112014098966119-pat00038
를 통해 도출할 수 있다.
즉, 상기한 [수학식 11]과 [수학식 12]를 [수학식 10]에 대입하고 [수학식 8]의
Figure 112014098966119-pat00039
를 이용하면, 이하의 [수학식 14]와 같이 하여
Figure 112014098966119-pat00040
의 1차 확장식인
Figure 112014098966119-pat00041
를 구할 수 있다.
마찬가지로,
Figure 112014098966119-pat00042
의 2차 확장식
Figure 112014098966119-pat00043
또한 상기한 [수학식 11]과 [수학식 12]를 [수학식 10]에 대입하고 [수학식 9]의
Figure 112014098966119-pat00044
를 이용하여, 이하의 [수학식 15]와 같이 계산할 수 있다.
[수학식 14]
Figure 112014098966119-pat00045

여기서, 나머지 성분들은 0(zero)이다.
[수학식 15]
Figure 112014098966119-pat00046

여기서, 나머지 성분들은 0(zero)이다.
Figure 112014098966119-pat00047
행렬의 대각 성분들은 P파와 S파의 전파에 관여하고, 나머지 성분들(off-diagonal entries)은 모드 결합(mode coupling)에 관여하게 되며, Fourier 영역에서 공간 변수 (ε)를 분리하지 않으면 모든 계산 노드에서 Fourier 변환을 실시하여 공간과 파수 영역을 오가야 하는 어려움이 있다.
그러나 EGS 전파자의 경우, 공간 변수를 효율적으로 Fourier 영역에서 분리하여 계산하기 때문에 노드의 수만큼 Fourier 변환을 하는 것이 아니라 확장된 항(term)의 수만큼만 변환을 실시하게 되어 파의 전파에 관한 계산을 빠르게 수행할 수 있다.
이와 같이 공간 변수를 주파수-파수 영역으로부터 효율적으로 분리하기 위해,
Figure 112014098966119-pat00048
는, 이하의 [수학식 16]에 나타낸 바와 같이, 위상 항(αv)과 공간변수 항(xμ)의 곱 형태로 전개할 수 있다(Le Rousseau, 2001, 참고문헌 9 참조).
[수학식 16]
Figure 112014098966119-pat00049

여기서, λ는 [수학식 14] 및 [수학식 15]에 있는 항(term)의 수를 나타낸다.
따라서 [수학식 11] 내지 [수학식 13] 및 [수학식 16]을 상기한 [수학식 1]과 [수학식 2]에 대입하면, 이하의 [수학식 17]과 같은 최종적인 EGS 전파자를 도출할 수 있다.
[수학식 17]
Figure 112014098966119-pat00050

여기서,
Figure 112014098966119-pat00051
는 스크린(screen) 항으로 주파수-공간(f-x) 영역에서 계산되고,
Figure 112014098966119-pat00052
는 위상 이동(phase-shift) 항으로 주파수-파수(f-k) 영역에서 계산이 수행된다.
또한, N은 정규화 연산자(normalizing operator)로 무한급수 전개시 절단오차로 인해 파수 또는 전파각에 따라 에너지가 부정확하게 증폭 또는 감쇠되는 현상을 방지하기 위해 사용된다(Le Rousseau and De Hoop, 2003, 참고문헌 8 참조).
즉, 도 1을 참조하면, 도 1은 EGS 전파자의 개념을 나타내는 도면으로, 상기한 [수학식 17]의 계산을 개념적으로 나타내는 도면이다.
즉, 도 1에 나타낸 바와 같이, P파와 S파로 분리된 파동장은 각각 두 개의 배경 매질(P파와 S파의 배경 속도가 각각
Figure 112014098966119-pat00053
인 매질과
Figure 112014098966119-pat00054
인 매질의 얇은 막)로 진행한다.
S파가 진행하는 두 번째
Figure 112014098966119-pat00055
매질의 경우, P파의 배경속도로
Figure 112014098966119-pat00056
가 아닌
Figure 112014098966119-pat00057
를 사용하는 것은 S파의 속도가 P파의 배경 속도보다 낮기 때문에 브랜치 포인트(branch point)가 발생할 가능성이 있기 때문이다.
즉, 도 1에서 있어서, P파는
Figure 112014098966119-pat00058
의 배경속도를 가진 매질을 진행하게 되는데 이때 발생하는 S파(그림에서
Figure 112014098966119-pat00059
)는 모드 전환에만 기여를 하고 실제 다음 얇은 막에서는 사용되지 않는다.
마찬가지로, S파는
Figure 112014098966119-pat00060
의 배경속도를 가진 매질을 전파하는데, 이때 발생한 P파(
Figure 112014098966119-pat00061
)는 모드 전환에만 사용이 되고 다음 막에서는 사용되지 않는다.
즉,
Figure 112014098966119-pat00062
의 배경속도를 가진 매질을 전파한 P파와
Figure 112014098966119-pat00063
의 배경속도를 가진 매질을 전파한 S파가 실제 다음 막으로 진행하는 P파와 S파가 된다.
또한, 상기한 바와 같이, 얇은 막을 통과한 P파와 S파 파동장은 M0에 의해 두 모드가 합쳐져 원래의 파동장(Vx, Vz)이 되고, 공간 영역(F-X 영역)에서
Figure 112014098966119-pat00064
에 의해 스크린 항(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 112014098966119-pat00065

여기서,
Figure 112014098966119-pat00066
는 푸리에 영역에서의 스칼라(scalar) 음원 파동장(source wavefield)을 나타내고,
Figure 112014098966119-pat00067
은 푸리에 영역에서의 스칼라 수진기 파동장(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 112014098966119-pat00068

여기서, 분모항에 나타나 있는 ε는 수진기 파동장의 절대값의 제곱에 일정 스케일러(scaling fraction)(0 < λ < 1)를 적용한 값으로
Figure 112014098966119-pat00069
로 표현할 수 있으며, 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 112014098966119-pat00070

여기서, γ는 반사지점에서의 반사각을 나타내며 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 112014098966119-pat00071
을 일괄적으로 곱하여 생성하였으며, 밀도 모델은 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 구조보정 방법.

    Figure 112015060550725-pat00090


    (여기서, Iij는 최종 영상화 이미지,
    Figure 112015060550725-pat00091
    는 푸리에 영역에서의 스칼라(scalar) 음원 파동장(source wavefield),
    Figure 112015060550725-pat00092
    은 푸리에 영역에서의 스칼라 수진기 파동장(receiver wavefield)을 나타내고, ε는
    Figure 112015060550725-pat00093
    로 표현할 수 있으며, i 와 j 는 음원과 수진기의 벡터 파동장(vector wavefield)으로, 다성분 탄성파 자료에서 수평성분(x-component)과 수직성분(z-component)을 각각 의미함.)
  2. 제 1항에 있어서,
    상기 EGS 전파자는,
    이하의 수학식에 나타낸 EGS 전파자를 이용하도록 구성되는 것을 특징으로 하는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법.

    Figure 112014098966119-pat00072


    (여기서, xμ(μ = 1, 2)과 x3는 각각 수평과 수직 좌표를 나타내고, s = -iω(ω는 각주파수) 이며, αv = kv/iω(kv는 수평성분 파수, v = 1, 2)이고,
    Figure 112014098966119-pat00073
    와
    Figure 112014098966119-pat00074
    는 각각 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. 삭제
  7. 제 1항에 있어서,
    상기 방법은,
    이하의 수학식을 이용하여 주파수-파수 영역에서 반사 지점에서의 반사각을 구하는 것에 의해, 매질 경계면의 반사점에서 극성이 변화하게 되는 S파 극성 역전 현상을 보정하는 극성보정단계를 더 포함하여 구성되는 것을 특징으로 하는 탄성파 다성분 자료에 대한 중합전 EGS 구조보정 방법.

    Figure 112015060550725-pat00079


    (여기서, γ는 반사지점에서의 반사각을 나타내며, kh와 kz는 각각 거리 방향의 파수와 심도 방향의 파수를 나타냄.)

  8. 삭제
KR1020140140533A 2014-10-17 2014-10-17 탄성파 다성분 자료에 대한 중합전 ｅｇｓ 구조보정 방법 KR101549388B1 (ko)

Priority Applications (3)

Application Number Priority Date Filing Date Title
KR1020140140533A KR101549388B1 (ko) 2014-10-17 2014-10-17 탄성파 다성분 자료에 대한 중합전 ｅｇｓ 구조보정 방법
PCT/KR2015/010951 WO2016060513A1 (ko) 2014-10-17 2015-10-16 탄성파 다성분 자료에 대한 중합전 egs 구조보정 방법
US15/508,957 US20170299745A1 (en) 2014-10-17 2015-10-16 Prestack egs migration method for seismic wave multi-component data

Applications Claiming Priority (1)

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

Publications (1)

Publication Number Publication Date
KR101549388B1 true KR101549388B1 (ko) 2015-09-02

Family

ID=54246978

Family Applications (1)

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

Country Status (3)

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

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105629301A (zh) * 2016-03-29 2016-06-01 中国地质大学（北京） 薄层弹性波反射系数快速求解方法
KR101865016B1 (ko) * 2018-01-16 2018-06-05 채휘영 지표투과레이다 자료 기반의 공동 검출 방법 및 공동 검출 장치
KR101907635B1 (ko) 2016-05-30 2018-10-12 한양대학교 산학협력단 역시간 구조보정을 이용하는 지하구조 영상장치 및 방법
CN110058300A (zh) * 2018-10-30 2019-07-26 南方科技大学 一次波重建方法、装置、终端设备及存储介质
CN113406698A (zh) * 2021-05-24 2021-09-17 中国石油大学(华东) 一种基于纵横波解耦的双相介质弹性波逆时偏移成像方法
CN113933900A (zh) * 2021-10-15 2022-01-14 张焕钧 一种隧道超前探测成像方法及装置
CN115993650A (zh) * 2023-03-22 2023-04-21 中国石油大学(华东) 一种基于棱柱波的地震干涉成像方法
CN117518265A (zh) * 2023-11-01 2024-02-06 大连理工大学 层状介质场地下基于频散特性的多模态面波波场反演方法

Families Citing this family (20)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105938203B (zh) * 2016-06-24 2018-07-10 中国石油天然气股份有限公司 一种储层特性的检测方法及装置
CN108072896B (zh) * 2016-11-18 2019-08-27 中国石油化工股份有限公司 一种全自动地震波初至拾取方法及系统
CN106772585B (zh) * 2017-01-26 2018-11-09 中国科学院地质与地球物理研究所 一种基于弹性波解耦方程的优化拟解析方法及装置
CN106896403B (zh) * 2017-05-05 2019-09-13 中国石油天然气集团有限公司 弹性高斯束偏移成像方法和系统
CN107918156B (zh) * 2017-10-30 2019-09-06 中国石油天然气集团公司 检测海底节点采集地震数据极性的方法及装置
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 中南大学 基于主频波数域均值滤波的地震数据去噪方法及系统
CN116413801B (zh) * 2023-02-13 2023-09-29 中国石油大学(北京) 一种各向异性介质弹性波高精度成像方法和系统
CN117233838B (zh) * 2023-09-20 2024-04-05 长江大学 一种二维vti介质中的弹性准纵横波场分离和逆时偏移成像方法
CN117249894B (zh) * 2023-11-16 2024-04-05 自然资源部第一海洋研究所 一种水下远场声传播在海底透射厚度的诊断方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR101355106B1 (ko) 2011-06-27 2014-01-22 서울대학교산학협력단 복소주파수 그룹을 기초로 하는 지하구조 영상화 방법

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7797110B2 (en) * 2007-06-26 2010-09-14 Shin's Geophysics Method for velocity analysis using waveform inversion in Laplace domain for geophysical imaging
US8537638B2 (en) * 2010-02-10 2013-09-17 Exxonmobil Upstream Research Company Methods for subsurface parameter estimation in full wavefield inversion and reverse-time migration
EP2691795A4 (en) * 2011-03-30 2015-12-09 CONVERGENCE SPEED OF COMPLETE WAVELENGTH INVERSION USING SPECTRAL SHAPING

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR101355106B1 (ko) 2011-06-27 2014-01-22 서울대학교산학협력단 복소주파수 그룹을 기초로 하는 지하구조 영상화 방법

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
논문 2012.08

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105629301A (zh) * 2016-03-29 2016-06-01 中国地质大学（北京） 薄层弹性波反射系数快速求解方法
CN105629301B (zh) * 2016-03-29 2018-02-09 中国地质大学（北京） 薄层弹性波反射系数快速求解方法
KR101907635B1 (ko) 2016-05-30 2018-10-12 한양대학교 산학협력단 역시간 구조보정을 이용하는 지하구조 영상장치 및 방법
KR101865016B1 (ko) * 2018-01-16 2018-06-05 채휘영 지표투과레이다 자료 기반의 공동 검출 방법 및 공동 검출 장치
CN110058300A (zh) * 2018-10-30 2019-07-26 南方科技大学 一次波重建方法、装置、终端设备及存储介质
CN113406698A (zh) * 2021-05-24 2021-09-17 中国石油大学(华东) 一种基于纵横波解耦的双相介质弹性波逆时偏移成像方法
CN113406698B (zh) * 2021-05-24 2023-04-21 中国石油大学(华东) 一种基于纵横波解耦的双相介质弹性波逆时偏移成像方法
CN113933900A (zh) * 2021-10-15 2022-01-14 张焕钧 一种隧道超前探测成像方法及装置
CN115993650A (zh) * 2023-03-22 2023-04-21 中国石油大学(华东) 一种基于棱柱波的地震干涉成像方法
CN115993650B (zh) * 2023-03-22 2023-06-06 中国石油大学(华东) 一种基于棱柱波的地震干涉成像方法
CN117518265A (zh) * 2023-11-01 2024-02-06 大连理工大学 层状介质场地下基于频散特性的多模态面波波场反演方法

Also Published As

Publication number Publication date
WO2016060513A1 (ko) 2016-04-21
US20170299745A1 (en) 2017-10-19

Similar Documents

Publication Publication Date Title
KR101549388B1 (ko) 탄성파 다성분 자료에 대한 중합전 ｅｇｓ 구조보정 방법
US11092707B2 (en) Determining a component of a wave field
RU2587498C2 (ru) Инверсия одновременных источников для данных сейсмоприемной косы с взаимнокорреляционной целевой функцией
Brossier et al. Which data residual norm for robust elastic frequency-domain full waveform inversion?
KR101548976B1 (ko) 지진 표면파들의 파형들을 사용하는 토양 특성들의 추정
Yao et al. Separation of migration and tomography modes of full‐waveform inversion in the plane wave domain
KR102021276B1 (ko) 진폭 보존을 갖는 fwi 모델 도메인 각도 스택들
US10935680B2 (en) Generating geophysical images using directional oriented wavefield imaging
Xu et al. 3D common image gathers from reverse time migration
CA2972033C (en) Multistage full wavefield inversion process that generates a multiple free data set
CA2808083C (en) Hybrid method for full waveform inversion using simultaneous and sequential source method
Cao et al. 3-D multiparameter full-waveform inversion for ocean-bottom seismic data using an efficient fluid–solid coupled spectral-element solver
KR101820850B1 (ko) 직접 파형 역산의 반복 적용을 이용한 탄성파 영상화 장치 및 방법
Sandberg et al. Full-wave-equation depth extrapolation for migration
Lamert et al. Full waveform inversion for advance exploration of ground properties in mechanized tunneling
Feng et al. Multiscale phase inversion for vertical transverse isotropic media
KR20170009609A (ko) 반복적 직접 파형 역산 및 완전 파형 역산을 이용한 탄성파 영상화 장치 및 방법
US10379245B2 (en) Method and system for efficient extrapolation of a combined source-and-receiver wavefield
Jia et al. Superwide-angle one-way wave propagator and its application in imaging steep salt flanks
Kazei et al. Amplitude-based DAS logging: Turning DAS VSP amplitudes into subsurface elastic properties
Kim et al. Prestack elastic generalized-screen migration for multicomponent data
Gao et al. Waveform tomography of two-dimensional three-component seismic data for HTI anisotropic media
Przebindowska Acoustic full waveform inversion of marine reflection seismic data
Huang et al. Hybrid local Born/Rytov Fourier migration method
EP3092520A1 (en) Determining a component of a wave field

Legal Events

Date Code Title Description
E701 Decision to grant or registration of patent right
GRNT Written decision to grant
FPAY Annual fee payment

Payment date: 20180625

Year of fee payment: 4

FPAY Annual fee payment

Payment date: 20190626

Year of fee payment: 5