KR101318994B1 - Method and apparatus of estimating underground structure using a plurality of weighting value - Google Patents

Method and apparatus of estimating underground structure using a plurality of weighting value Download PDF

Info

Publication number
KR101318994B1
KR101318994B1 KR1020120019214A KR20120019214A KR101318994B1 KR 101318994 B1 KR101318994 B1 KR 101318994B1 KR 1020120019214 A KR1020120019214 A KR 1020120019214A KR 20120019214 A KR20120019214 A KR 20120019214A KR 101318994 B1 KR101318994 B1 KR 101318994B1
Authority
KR
South Korea
Prior art keywords
data
underground medium
steep slope
maximum
underground
Prior art date
Application number
KR1020120019214A
Other languages
Korean (ko)
Other versions
KR20130097504A (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 KR1020120019214A priority Critical patent/KR101318994B1/en
Publication of KR20130097504A publication Critical patent/KR20130097504A/en
Application granted granted Critical
Publication of KR101318994B1 publication Critical patent/KR101318994B1/en

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. analysis, for interpretation, for correction
    • G01V1/30Analysis
    • G01V1/301Analysis for determining seismic cross-sections or geostructures
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. analysis, for interpretation, for correction
    • G01V1/30Analysis
    • G01V1/306Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. analysis, for interpretation, for correction
    • G01V1/32Transforming one recording into another or one representation into another

Abstract

탄성파(인공지진파)를 이용하여 지하매질의 구조를 추정하는 기술이 개시된다. 본 발명에 따른 지하 매질을 통과하여 지표면에서 기록된 탄성파 데이터와 상기 가정한 초기 지하 매질에 대한 모델링 데이터를 가지고 파형역산 알고리즘을 적용하여 상기 지하 매질의 구조를 추정하는 방법은, 상기 측정된 시간 영역의 데이터를 소정의 변환 영역의 데이터로 변환하고, 상기 변환 영역에서의 상기 지하 매질 구조에 대한 측정 데이터와 모델링 데이터의 차이를 이용하여 각 변환변수에서의 목적 함수(objective function)를 정의하고 잔차를 계산하는 단계; 상기 각 변환변수에서의 목적 함수를 최소로 하는 최대 급경사 방향을 합산하여 전체 변환변수의 범위에 대한 최대 급경사 방향을 계산하되, 제1가중치를 적용하여 상기 각 변환변수에서의 최대 급경사 방향을 정규화(normalization)하고, 상기 정규화된 각 변환변수별 최대 급경사 방향에 제2가중치를 추가로 적용하여 합산하는 단계를 포함한다.A technique for estimating the structure of an underground medium using a seismic wave (artificial earthquake wave) is disclosed. The method for estimating the structure of the underground medium by applying a waveform inversion algorithm with the seismic data recorded on the surface of the ground and passing through the underground medium according to the present invention and modeling data for the assumed initial underground medium may include: Transforms the data of the transformed data into a data of a predetermined transformation region, and defines an objective function in each transformation variable by using the difference between the measured data and the modeling data of the structure of the underground medium in the transformation region. Calculating; Compute the maximum steep slope for the range of all the transformation variables by summing the maximum steep slopes that minimize the objective function in each of the transformation variables, and normalize the maximum steep slope in each transformation variable by applying a first weight value. normalizing) and additionally applying a second weighting value to the maximum steep slope for each of the normalized conversion variables.

Description

복수의 가중치를 이용한 지하매질 구조 추정 방법 및 장치 {Method and apparatus of estimating underground structure using a plurality of weighting value}Method and apparatus for estimating underground structure using a plurality of weights {a method and apparatus of estimating underground structure using a plurality of weighting value}

탄성파(인공지진파)를 이용하여 지하매질의 구조를 추정하는 기술과 관련된다.
Related to techniques for estimating the structure of underground media using seismic waves (artificial earthquake waves).

지하에 매장된 광물이나 원유 등 지하자원을 개발하는데는 많은 비용과 노력이 필요하다. 보다 구체적으로는 직간접적인 여러가지 방법의 지질탐사를 통해 지하자원이 매장되어 있을 가능성이 가장 높은 지역을 찾아내고, 그 중 몇 몇 지점에 대해 직접 시추 또는 간접적인 자료분석을 통해 지하자원의 유무와 매장량 등을 확인해 보는 과정을 거친다.Developing underground resources such as minerals and crude oil buried underground requires a lot of cost and effort. More specifically, the area most likely to be buried underground is found through direct and indirect geological exploration, and the presence and reserves of underground resources are obtained through direct drilling or indirect data analysis of some of them. Check the back.

한편, 지하자원의 유무, 위치 등은 지진파와 같은 탄성파를 사용하여 알아낼 수 있다. 즉, 인공적으로 지진파를 발생시키고, 지하매질을 통과하여 수신된 지진파 데이터를 분석함으로써 지하매질의 밀도 등을 알아낼 수 있고, 이를 통해 지하매질의 구조를 파악함으로써 지하자원의 유무, 매장 위치 등을 알아낼 수 있다. 이 과정에서 측정된 데이터와의 오차를 최소화시키는 모델링 데이터를 제공하는 지하모델을 구하는 일반적인 파형역산 알고리즘이 적용된다.On the other hand, the presence, location, etc. of the underground resources can be found using the seismic waves, such as seismic waves. In other words, by artificially generating seismic waves and analyzing the seismic wave data received through the underground media, the density of the underground media can be determined, and the structure of the underground media can be used to determine the presence or absence of underground resources and burial location. Can be. In this process, a general waveform inversion algorithm is applied to obtain an underground model that provides modeling data that minimizes errors from measured data.

일반적으로는 실제 측정된 탄성파 데이터와, 현장 지질답사나 시추 코어 분석등을 통해 얻어진 사전 정보를 종합해 가정한 초기 지하매질(초기모델)에 대한 모델링 데이터를 가지고 파형역산을 수행하는데, 파형역산에 필요한 여러 파라미터를 반복적으로 갱신해가면서 지하매질에 대한 지진파의 속도와 지하매질의 밀도를 함께 계산한다. In general, waveform inversion is performed with modeling data for the initial underground medium (initial model), which is a combination of actual measured seismic data and preliminary information obtained through field geological surveys or drilling core analysis. By repeatedly updating the necessary parameters, the velocity of the seismic wave with respect to the underground medium and the density of the underground medium are calculated together.

그러나, 파형역산 알고리즘은 아직 많은 제한점을 가지고 있으며, 그 중 하나가 역산 결과가 초기 모델에 매우 민감하다는 것이다. 특히 대표적인 석유 유망구조인 암염 모델의 경우에, 초기 모델이 좋지 못할 경우 역산 결과 추정된 지하매질 구조의 정확도가 떨어지는데, 제한된 사전 정보만으로 좋은 초기 모델을 정하는 것은 매우 어렵다. 따라서 역산 결과 추정된 지하매질 구조의 정확도를 향상시키기 위해 보다 나은 초기모델을 정하기 위한 연구가 진행되어 왔다.
However, waveform inversion algorithms still have many limitations, one of which is that inversion results are very sensitive to the initial model. Particularly in the case of the rock salt model, which is a representative petroleum promising structure, if the initial model is not good, the accuracy of the underground medium structure estimated by the inversion is inferior. Therefore, studies have been conducted to establish a better initial model to improve the accuracy of the underground media structure estimated by the inversion.

목적함수를 최소로 하는 최대 급경사 방향(gradient) 계산시, 각 변환변수별로 계산된 최대 급경사 방향은 서로 다른 공간적 분해능(spatial resolution)을 가지는데 보다 정확하고 효율적으로 최대 급경사 방향을 결정하기 위해 이들 각각은 송신원이 제거된 잔여파동장(residual wavefield)의 특성을 반영하여 가중되어야 한다. 따라서 각 변환변수별 최대 급경사 방향을 합산하여 전체 변환변수의 범위에 대한 최대 급경사 방향을 계산하되, 제1가중치를 적용하여 각 변환변수에서의 최대 급경사 방향을 정규화(normalization)하고, 정규화된 각 변환변수별 최대 급경사 방향에 제2가중치를 더 적용하여, 보다 정확하게 지하매질 구조를 추정하는 방법 및 장치를 제공하고자 한다.
When calculating the maximum gradient which minimizes the objective function, the maximum steep slope calculated for each transform variable has different spatial resolutions, each of which is used to determine the maximum steep slope more accurately and efficiently. Should be weighted to reflect the characteristics of the residual wavefield from which the source has been removed. Therefore, the maximum steep slope is calculated by summing the maximum steep slope directions of each transform variable, normalizing the maximum steep slope direction in each transform variable by applying a first weight value, and normalizing each normalized transform. The present invention provides a method and apparatus for more accurately estimating the structure of the underground medium by further applying a second weighting value to the maximum steep slope of each variable.

본 발명의 일 양상에 따른, 실제 지하 매질을 통과하여 지표면에서 기록된 탄성파 데이터와 상기 가정한 초기 지하 지질구조 모델에 대한 모델링 데이터를 가지고 파형역산 알고리즘을 적용하여 상기 실제 지하 매질의 구조를 추정하는 방법은, 상기 측정된 시간 영역의 데이터를 소정의 변환 영역의 데이터로 변환하고, 상기 변환 영역에서의 상기 지하 매질 구조에 대한 측정 데이터와 모델링 데이터의 차이를 이용하여 각 변환변수에서의 목적 함수(objective function)를 정의하고 잔차를 계산하는 단계; 상기 각 변환변수에서의 목적 함수를 최소로 하는 최대 급경사 방향을 합산하여 전체 변환변수의 범위에 대한 최대 급경사 방향을 계산하되, 제1가중치를 적용하여 상기 각 변환변수에서의 최대 급경사 방향을 정규화(normalization)하고, 상기 정규화 된 각 변환변수별 최대 급경사 방향에 제2가중치를 추가로 적용하여 합산하는 단계를 포함한다.According to an aspect of the present invention, the structure of the actual underground medium is estimated by applying a waveform inversion algorithm with the seismic data recorded on the surface of the earth through the actual underground medium and the modeling data for the assumed initial underground geologic model. The method converts the measured time domain data into data of a predetermined transformation region and uses the difference between the measured data and the modeling data for the underground medium structure in the transformation region to obtain an objective function ( defining an objective function) and calculating a residual; Compute the maximum steep slope for the range of all the transformation variables by summing the maximum steep slopes that minimize the objective function in each of the transformation variables, and normalize the maximum steep slope in each transformation variable by applying a first weight value. and normalizing the second weighted value by additionally applying a second weight value to a maximum steep slope for each normalized conversion variable.

이때 상기 변환 영역은 주파수 도메인이고, 상기 제1가중치는 상기 지하 매질 구조를 나타내는 소정의 파라미터(pω)에 대한 상기 각 변환변수에서 계산된 최대 급경사 방향 벡터에서 그 절대값이 가장 큰 값의 역수로 정의될 수 있다.In this case, the transform region is a frequency domain, and the first weight value is an inverse of a value whose largest absolute value is in the maximum steepness direction vector calculated from each transform variable for a predetermined parameter p ω representing the underground medium structure. It can be defined as.

또한 상기 변환 영역은 주파수 도메인이고, 상기 제2가중치는 상기 탄성파 송신원이 제거된 측정 데이터와 모델링 데이터 사이의 잔여파동장을 송신원으로 하여 계산된 역전파된 파동장의 진폭의 평균으로 정의될 수 있다.In addition, the transform region is a frequency domain, and the second weighting value may be defined as an average of amplitudes of the back propagated wave fields calculated using a residual wave field between the measured data from which the acoustic wave transmitter is removed and modeling data as a transmission source.

또한 본 발명의 다른 양상에 따른 실제 지하 매질을 통과하여 지표면에서 기록된 탄성파 데이터와 상기 가정한 초기 지하 지질구조에 대한 모델링 데이터를 가지고 파형역산 알고리즘을 적용하여 상기 실제 지하 매질의 구조를 추정하는 장치는, 상기 측정된 시간 영역의 데이터를 소정의 변환 영역의 데이터로 변환하는 변환부; 상기 추정하고자 하는 지하 매질을 모델링하기 위해 초기 모델을 설정하고, 상기 변환 영역에서의 모델링 데이터를 산출하는 모델링 데이터 생성부; 상기 각 변환 변수에서의 목적 함수를 정의하고 측정 데이터와 모델링 데이터 사이의 잔차를 계산하는 잔차 산출부; 상기 각 변환변수에서의 최대 급경사 방향을 독립적으로 계산한 뒤 합산하여 전체 변환변수의 범위에 대한 최대 급경사 방향을 계산하되, 제1가중치를 적용하여 상기 각 변환변수에서의 최대 급경사 방향을 정규화(normalization)하고, 상기 정규화된 각 변환변수별 최대 급경사 방향에 제2가중치를 더 적용하여 합산하여 최대 급경사 방향을 계산하는 최대 급경사 방향 산출부; 및 상기 전체 변환 변수 범위에서 계산된 최대 급경사 방향을 상기 목적 함수가 최소가 되는 방향으로 상기 모델 파라미터를 업데이트하는 갱신부를 포함한다.
In addition, the apparatus for estimating the structure of the actual underground medium by applying a waveform inversion algorithm with the seismic data recorded on the surface of the earth surface and modeling data for the assumed initial underground geological structure according to another aspect of the present invention The conversion unit may include: a conversion unit converting the measured time domain data into data of a predetermined conversion region; A modeling data generator configured to set an initial model to model the underground medium to be estimated, and to calculate modeling data in the transform region; A residual calculator for defining an objective function in each of the conversion variables and calculating a residual between the measured data and the modeling data; Calculate the maximum steep slope for the entire range of the transformed variable by independently calculating and summing up the maximum steep slope in each of the transform variables, and normalizing the maximum steep slope in each of the transform variables by applying a first weight value. A maximum steep direction calculation unit configured to calculate a maximum steep slope direction by adding a second weighting value to the normalized maximum steep slope direction for each of the normalized conversion variables; And an updating unit for updating the model parameter in a direction in which the objective function is minimum in the maximum steep slope direction calculated in the entire conversion variable range.

제1가중치와 제2가중치를 함께 적용함으로써, 초기 지하모델이 부정확하게 추정된 경우에도 보다 정확하게 지하매질 구조를 추정할 수 있다.
By applying the first and second weights together, it is possible to estimate the structure of the underground medium more accurately even when the initial underground model is incorrectly estimated.

도 1은 모델링을 위해 사용하는 송신원의 진폭을 주파수 영역에서 도시한 도면이다.
도 2는 본 발명의 지하 매질 구조 추정방법을 테스트하기 위한 단순화된 모델을 도시한 도면이다.
도 3은 시간 영역에서 송신원이 포함된 잔여파동장을 구한 도면이다.
도 4는 각 주파수별로 잔차의 제곱으로 정의한 오차를 출력한 도면이다.
도 5는 본 발명의 일실시예에 따른 지하매질 구조 추정 방법의 흐름도이다.
도 6은 본 발명의 일실시예에 따른 지하매질 구조 추정 장치의 구성도이다.
도 7은 방법1을 이용해 각각의 모델에 대한 탄성 파형역산을 수행한 결과, 라메상수(λ)에 대해 첫 번째 이터레이션에서 계산된 최대 급경사 방향을 도시한 도면이다.
도 8은 방법2를 이용해 각각의 모델에 대한 탄성 파형역산을 수행한 결과, 라메상수(λ)에 대해 첫 번째 이터레이션에서 계산된 최대급경사 방향을 도시한 도면이다.
도 9는 본 발명의 가중치기법을 이용해 각각의 모델에 대한 탄성 파형역산을 수행한 결과, 라메상수(λ)에 대해 첫 번째 이터레이션에서 계산된 최대 급경사 방향을 도시한 도면이다.
도 10은 방법1을 이용해 파형역산을 수행하여 추정한 지하매질에서의 P파 속도를 도시한 도면이다.
도 11은 방법2를 이용해 파형역산을 수행하여 추정한 지하매질에서의 P파 속도를 도시한 도면이다.
도 12는 가중치기법을 이용해 파형역산을 수행하여 추정한 지하매질에서의 P 속도를 도시한 도면이다.
도 13은 가정한 모델링 영역의 정 중앙(3km)에서 획득한 P파 속도 단면을 도시한 도면이다.
도 14는 파형역산이 끝난 후, 최종적으로 송신원이 제거된 RMS(Root Mean Square) 오차를 주파수 별로 출력한 도면이다.
도 15는 SEG/EAGE 암염모델의 P파 속도와 S파 속도 모델로 탄성 파형역산에서 실제 지질모델로 가정한 모델을 출력한 도면이다.
도 16은 초기모델의 P파 속도가 1.5 ~ 3.06 km로 선형적으로 증가하도록 설정했을 때 방법1을 이용해 암염모델에 대한 탄성 파형역산을 수행한 결과 추정한 P파 속도와 S파 속도를 도시한 도면이다.
도 17은 방법1을 이용해 암염모델에 대한 탄성 파형역산을 수행한 결과 파라미터 라메상수(λ와 μ)에 대해 300번째 이터레이션에서 계산된 최대 급경사 방향을 도시한 도면이다.
도 18은 초기모델의 P파 속도가 1.5 ~ 3.06 km로 선형적으로 증가하도록 설정했을 때 가중치기법을 이용해 암염모델에 대한 탄성 파형역산을 수행한 결과 추정한 P파 속도와 S파 속도를 도시한 도면이다.
도 19는 가중치기법을 이용해 암염모델에 대한 탄성 파형역산을 수행한 결과 파라미터 라메상수(λ와 μ)에 대해 300번째 이터레이션에서 계산된 최대 급경사 방향을 도시한 도면이다.
도 20은 각각의 방법을 이용해 암염모델에 대한 탄성 파형역산을 수행한 결과 추정한 6 km와 9 km 지점에서의 P파 속도 단면을 도시한 도면이다.
1 is a diagram showing the amplitude of a transmission source used for modeling in the frequency domain.
2 illustrates a simplified model for testing the method for estimating the underground medium structure of the present invention.
3 is a diagram of a residual wave field including a transmission source in a time domain.
4 is a diagram showing an error defined by the square of the residual for each frequency.
5 is a flowchart of a method for estimating an underground medium structure according to an embodiment of the present invention.
Figure 6 is a block diagram of a structure for estimating the underground medium according to an embodiment of the present invention.
FIG. 7 is a diagram showing the maximum steep slope calculated in the first iteration for the lame constant λ as a result of performing the elastic waveform inversion for each model using Method 1. FIG.
FIG. 8 is a diagram showing the maximum steep slope calculated in the first iteration for the lame constant λ as a result of performing the elastic waveform inversion of each model using the method 2. FIG.
FIG. 9 is a diagram showing the maximum steep slope calculated in the first iteration for the lame constant λ as a result of performing the elastic waveform inversion of each model using the weighting method of the present invention.
FIG. 10 is a diagram showing P wave velocity in an underground medium estimated by performing waveform inversion using Method 1. FIG.
FIG. 11 is a diagram showing P wave velocity in an underground medium estimated by performing waveform inversion using Method 2. FIG.
FIG. 12 is a diagram illustrating P velocity in an underground medium estimated by performing waveform inversion using a weighting technique.
FIG. 13 is a diagram showing a P-wave velocity cross section obtained at the center (3 km) of the assumed modeling region.
FIG. 14 is a diagram illustrating output of a root mean square (RMS) error for each frequency, from which a transmission source is finally removed after waveform inversion is completed.
FIG. 15 is a view showing a model assumed as a real geological model in elastic waveform inversion as a P-wave velocity and an S-wave velocity model of the SEG / EAGE rock salt model.
FIG. 16 shows the P-wave velocity and the S-wave velocity estimated by performing the elastic waveform inversion of the rock salt model using Method 1 when the P-wave velocity of the initial model was set to increase linearly from 1.5 to 3.06 km. Drawing.
FIG. 17 is a diagram showing the maximum steep slope calculated at the 300 th iteration for the parameter lame constants (λ and μ) as a result of performing elastic waveform inversion on the rock salt model using Method 1. FIG.
FIG. 18 illustrates the P-wave velocity and the S-wave velocity estimated by performing the elastic waveform inversion of the rock salt model using the weighting method when the P-wave velocity of the initial model is set to increase linearly from 1.5 to 3.06 km. Drawing.
FIG. 19 is a diagram showing the maximum steep slope calculated at the 300th iteration for the parameter lame constants (λ and μ) as a result of performing the elastic waveform inversion on the rock salt model using the weighting technique.
FIG. 20 is a diagram showing P-wave velocity cross-sections estimated at 6 km and 9 km as a result of performing elastic waveform inversion on the rock salt model using the respective methods.

이하 첨부된 도면을 참조하여 본 발명의 바람직한 실시예에 대해 상세히 설명한다. 본 발명을 설명함에 있어 관련된 공지 기능 또는 구성에 대한 구체적인 설명이 본 발명의 요지를 불필요하게 흐릴 수 있다고 판단되는 경우에는 그 상세한 설명을 생략할 것이다. 또한, 후술되는 용어들은 본 발명에서의 기능을 고려하여 정의된 용어들로서 이는 사용자, 운용자의 의도 또는 관례 등에 따라 달라질 수 있다. 그러므로 그 정의는 본 명세서 전반에 걸친 내용을 토대로 내려져야 할 것이다.DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS Reference will now be made in detail to the preferred embodiments of the present invention, examples of which are illustrated in the accompanying drawings. In the following description of the present invention, a detailed description of known functions and configurations incorporated herein will be omitted when it may make the subject matter of the present invention rather unclear. In addition, the terms described below are defined in consideration of the functions of the present invention, and this may vary depending on the intention of the user, the operator, or the like. Therefore, the definition should be based on the contents throughout this specification.

본 발명의 실시예에 따른 '파형 역산(waveform inversion)'이란, 현장에서 실제 측정된 데이터를 이용하여 특정지역의 지하 구조에 관한 정보(예컨대, 측정 대상 지역에 대한 속도 모델 또는 밀도 모델)를 유추하는 과정을 말한다. 이러한 파형 역산은 해석자가 임의의 지하구조 모델을 설정한 후, 설정된 모델에 대한 이론값을 구하는 모델링(modeling) 과정을 수반할 수 있다.The term 'waveform inversion' according to an embodiment of the present invention is to infer information (for example, a velocity model or a density model of a region to be measured) about the underground structure of a specific region using data measured in the field. The process of doing that. Such waveform inversion may involve a modeling process in which an analyst sets an arbitrary underground model and then obtains theoretical values for the set model.

예컨대, 본 발명의 실시예에 따른 파형 역산을 이용하여 지하 구조를 영상화 하는 경우, 모델링을 거쳐 계산된 이론값(이론 데이터)들과 실제 현장 탐사를 통해 얻어진 측정값(측정 데이터)들을 비교하여 얻어지는 차이값을 이용하여 새로운 지하 구조 모델을 만든다. 그리고 새로운 지하 구조 모델에 대한 모델링을 통해 구한 이론값들(모델링 데이터)을 다시 측정값들과 비교하는 과정을 반복적으로 수행한다. 이 경우에 이론값과 측정값의 비교 및 이를 통한 지하 구조 모델의 업데이트 과정은 그 차이값 또는 오차가 최소가 되거나 또는 소정의 임계치 이하가 될 때까지 반복될 수 있다. 이러한 차이값 또는 오차가 미리 결정된 소정의 범위 이내가 되면, 최종적으로 실제 지하 구조와 동일하거나 또는 유사한 지하 구조 모델을 얻을 수 있다.For example, when imaging the underground structure using the waveform inversion according to an embodiment of the present invention, the theoretical values calculated through modeling are compared with the measured values (measurement data) obtained through field surveys. Create a new underground structural model using the differences. And the process of comparing the theoretical values (modeling data) obtained through the modeling of the new underground structural model with the measured values is repeated. In this case, the comparison between the theoretical and measured values and the update of the underground structural model through the same may be repeated until the difference or error becomes minimum or below a predetermined threshold. When this difference or error falls within a predetermined range, it is possible to finally obtain an underground structural model that is the same as or similar to the actual underground structure.

이러한 본 발명의 실시예에 따른 파형 역산은, 측정 대상 지역의 지하 구조를 영상화하기 위한 영상 데이터를 생성하기 위해 각종 신호를 처리하는 계산 장치, 신호 처리 알고리즘이 기록된 컴퓨터로 읽을 수 있는 기록 매체, 이러한 계산 장치 또는 기록 매체 등을 통해 지하 구조를 영상화하는 방법 등에 의해 구체화될 수 있다.The waveform inversion according to the embodiment of the present invention includes a computing device for processing various signals to generate image data for imaging the underground structure of the measurement target area, a computer-readable recording medium having recorded thereon a signal processing algorithm, It can be embodied by a method of imaging an underground structure through such a computing device or a recording medium.

이하에서서는 본 발명에서 사용되는 파형역산, 가중치 기법에 대해 상세히 설명한다.Hereinafter, the waveform inversion and weighting technique used in the present invention will be described in detail.

(1)(One) 역산 이론Inversion theory

탄성파 파형 역산은 탄성파 탐사를 통해 기록한 자료와 동일한 모델링 데이터를 만들어내는 지질모델을 찾음으로써 실제 지하모델을 추정하는 방법으로 다음과 같은 수학적인 과정을 거쳐서 수행된다.Seismic waveform inversion is a method of estimating the actual underground model by finding a geological model that produces the same modeling data as the data recorded through seismic exploration.

실제 측정 데이터(d s )와 모델링 데이터(u s ) 사이의 차이를 정의하기 위해 l 2 norm 을 사용하면, 목적함수(오차함수)는 다음과 같이 정의된다. Using l 2 norm to define the difference between the actual measured data ( d s ) and the modeling data ( u s ), the objective function (error function) is defined as

Figure 112012015369898-pat00001
(1)
Figure 112012015369898-pat00001
(One)

이 때, s는 송신원, ω는 주파수, p는 추정하고자 하는 파라미터를 의미한다. 목적함수를 최소로 만드는 파라미터의 최대 급경사 방향을 구하기 위해 위의 식을 파라미터에 대해서 미분하면, 최대 급경사 방향은 다음과 같이 정의된다.In this case, s denotes a transmission source, ω denotes a frequency, and p denotes a parameter to be estimated. Differentiating the above equations for the parameters to find the maximum steepness of the parameter that minimizes the objective function, the maximum steepness is defined as

Figure 112012015369898-pat00002
(2)
Figure 112012015369898-pat00002
(2)

이 때, J k는 모델링을 위해 격자화된 각각의 셀(또는 격자)에서의 각 수진기에 기록된 데이터에 대한 민감도를 의미하는 자코비안 행렬이다. 자코비안 행렬을 직접 계산하는 것은 오랜 시간과 컴퓨터 메모리가 요구되므로, 자코비안 행렬을 계산하지 않고 최대 급경사 방향을 효율적으로 계산하기 위해 역전파 알고리즘이 널리 사용된다. 역전파 알고리즘을 이용하면, 최대 급경사 방향은 다음과 같이 계산할 수 있다.Where J k is a Jacobian matrix representing the sensitivity to the data recorded in each receiver in each cell (or grating) that has been truncated for modeling. Since direct computation of Jacobian matrices takes a long time and computer memory, backpropagation algorithms are widely used to efficiently calculate the maximum steep slope without computing the Jacobian matrices. Using the backpropagation algorithm, the maximum steep slope can be calculated as follows.

Figure 112012015369898-pat00003
(3)
Figure 112012015369898-pat00003
(3)

fv s,k는 가상 송신원 벡터를 의미하며 S는 임피던스 행렬이다. 이러한 방법을 그래디언트법(gradient method) 또는 최대 급경사법(steepest-descent method)이라고 부르는데, 탄성파 탐사에서 송신원이 표면에 위치하기 때문에 하부의 최대 급경사 방향(gradient)은 값이 매우 작다. 따라서 그래디언트법으로는 하부의 구조를 영상화하기에 어려움이 있다. 상부에 집중되어 있는 최대 급경사 방향값의 크기를 작게 하고, 하부의 최대 급경사 방향값을 키우기 위해서 헤시안(Hessian) 행렬를 이용해 최대 급경사 방향을 스케일링(scaling)하는 방법이 가우스-뉴튼법(Gauss-Newton method)이며, 가우스-뉴튼법을 이용해 계산한 최대 급경사 방향은 다음 식과 같이 표현된다.f v s, k denotes a virtual source vector and S is an impedance matrix. This method is called the gradient method or the steepest-descent method. In the seismic survey, the maximum steep gradient at the bottom is very small because the source is located on the surface. Therefore, it is difficult to image the structure of the lower part by the gradient method. In order to reduce the magnitude of the maximum steep slope value concentrated in the upper part and to increase the maximum steep slope value in the lower part, a method of scaling the maximum steep slope direction using a Hessian matrix is a Gauss-Newton method. method), and the maximum steep slope calculated using the Gauss-Newton method is expressed by the following equation.

Figure 112012015369898-pat00004
(4)
Figure 112012015369898-pat00004
(4)

하지만 자코비안 행렬을 계산하는 것은 큰 부담이 되고, 헤시안 행렬의 역행렬을 계산하는 것 또한 추가적인 컴퓨터 메모리와 계산시간이 너무 많이 소모되므로, 최대 급경사 방향 벡터를 스케일링하기 위해 가상송신원만으로 구성한 슈도헤시안(pseudo-Hessian) 행렬의 대각성분을 많이 사용하고 있으며 이는 다음 식과 같이 표현된다. However, computing Jacobian matrices is a heavy burden, and calculating the inverse of the Hessian matrix also takes too much additional computer memory and computation time, so pseudohessian consists of virtual sources only to scale the steepest gradient direction vector. The diagonal component of the (pseudo-Hessian) matrix is used a lot and is expressed as the following equation.

Figure 112012015369898-pat00005
(5)
Figure 112012015369898-pat00005
(5)

이 때, β는 헤시안 행렬의 스케일링 효과를 조절하기 위한 값이며 이 식을 통해 그래디언트 방향을 계산하면 모든 송신원에 대한 모든 주파수 성분의 목적함수에 대해서 오차를 최소화할 수 있는 파라미터 업데이트를 결정할 수 있다.
In this case, β is a value for adjusting the scaling effect of the Hessian matrix. By calculating the gradient direction through this equation, it is possible to determine the parameter update that can minimize the error for the objective function of all frequency components for all sources. .

(2) 일반적인 두 가지 스케일링 방법과 그 문제점 비교(2) Comparison of two common scaling methods and their problems

슈도헤시안 행렬을 이용한 스케일링은 두 가지 방법으로 수행할 수 있다. 첫 번째는 식 (6)과 같이 주파수루프 안에서 슈도헤시안 행렬을 적용하는 방법이며(이하 방법1이라고 한다), 두 번째는 식 (7)와 같이 주파수루프 밖에서 슈도헤시안 행렬을 적용하는 방법(이하 방법2라고 한다)이다.Scaling using a pseudohessian matrix can be performed in two ways. The first method is to apply the pseudohessian matrix in the frequency loop as in Eq. (6), and the second is to apply the pseudohesian matrix outside the frequency loop as in Equation (7). Hereinafter, method 2).

Figure 112012015369898-pat00006
(6)
Figure 112012015369898-pat00006
(6)

Figure 112012015369898-pat00007
(7)
Figure 112012015369898-pat00007
(7)

앞의 역산 이론에서 유도했던 것처럼 모든 송신원에 대한 모든 주파수 성분의 목적함수를 최소로 만드는 최대 급경사 방향을 구하는 것은 식 (7)에 따르는 것이다. 식 (6)은 각 주파수 성분에 대해서, 그 주파수 성분의 shot gather(공통 송신원 모음)에 대한 목적함수만을 최소로 만드는 최대 급경사 방향을 찾고, 이들을 모두 더함으로써, 최종적으로 하나의 최대 급경사 방향을 결정하게 된다. 이는 모든 주파수 성분에 대한 역산 이론과는 다른, 이론적으로 맞지 않는 방법이며, 이렇게 구한 최대 급경사 방향이 모든 주파수에 대한 오차를 최소화시킨다는 것을 보장할 수 없다. 그럼에도 불구하고 기존 연구에서는 방법1이 더 좋은 결과를 준다고 알려져 있다.As derived from the inversion theory above, finding the maximum steepest slope direction that minimizes the objective function of all frequency components for all sources is in accordance with Eq. (7). Equation (6) finds the maximum steepest slope for each frequency component by finding the maximum steepest direction that minimizes only the objective function for the shot gather of the frequency components, and adding them all together. Done. This is a theoretically inconsistent method, unlike the inversion theory for all frequency components, and we cannot guarantee that the maximum steepest direction thus obtained minimizes the error for all frequencies. Nevertheless, previous studies have found that Method 1 yields better results.

이는 방법2가 다음과 같은 문제점을 가지기 때문이다. 탄성파 탐사는 지표에서 다이너마이트와 같은 인위적인 송신원을 사용하여 신호를 발생시키고, 지하매질을 거쳐서 다시 표면으로 돌아온 파동장을 지표에서 기록하여 자료를 획득하며 이 자료에는 지하 매질의 물성에 대한 정보가 포함되어 있다. 탄성파 탐사에서 사용하는 송신원은 특정 주파수 성분이 지배적인 형태로, 모델링을 위해 사용하는 송신원은 주파수영역에서 도 1에 도시한 바와 같은 진폭을 가진다.This is because the method 2 has the following problems. Seismic exploration generates signals using an artificial source, such as dynamite, at the surface, and acquires data from the surface by recording wave fields that have returned to the surface through the underground medium, which includes information about the properties of the underground medium. have. The source used in the seismic survey is dominated by a specific frequency component, and the source used for modeling has an amplitude as shown in FIG. 1 in the frequency domain.

도 1은 모델링을 위해 사용하는 송신원의 주파수 영역에서의 파형을 도시한 도면이다.1 is a diagram illustrating waveforms in a frequency domain of a transmission source used for modeling.

가로축이 주파수, 세로축이 송신원의 진폭이며, 최대주파수(10 Hz)의 1/4 정도에서 최대 진폭을 가진다. 따라서 각 주파수 별로 계산되는 최대 급경사 방향도 최대주파수의 1/4 정도에서 가장 크게 계산되며, 매우 저주파 - 고주파인 자료는 값이 매우 작아서 최대 급경사 방향을 구성하는데 큰 기여를 하지 못한다. 따라서 송신원 효과를 제거해주는 기술을 쓰지 않는다면, 성공적인 역산을 기대할 수 없다. The horizontal axis represents the frequency and the vertical axis represents the amplitude of the transmission source, and has a maximum amplitude at about 1/4 of the maximum frequency (10 Hz). Therefore, the maximum steep slope calculated for each frequency is also calculated the most at about 1/4 of the maximum frequency, and very low-high frequency data is so small that it does not contribute much to construct the maximum steep slope. Therefore, unless you use a technique that eliminates the source effect, you cannot expect a successful inversion.

반면에 방법1에서는 슈도헤시안 행렬이 주파수루프 안에서 적용되며, 주파수 영역에서 송신원은 각 주파수에서 하나의 상수값을 가진다. 식 (6)에서 만약 송신원의 파형을 정확히 안다고 가정하면, 송신원이 가상송신원(F v s )에 한 번, 잔여파동장(u s -d s )에 한 번 곱해져 있기 때문에, 최종적으로 분자, 분모에 각각 두 번씩 곱해지게 되며 서로 상쇄되어 효율적으로 제거된다. 이는 송신원 파형역산이 성공적으로 수행되었을 때, 방법1을 사용하면 송신원의 효과를 제거하여 매우 저주파 ~ 고주파 성분의 자료가 역산에 기여할 수 있게 한다. 따라서 기존 연구에서 방법1이 송신원의 효과가 제거되지 않은 방법2보다 더 좋은 결과를 제공하였다. On the other hand, in method 1, the pseudohessian matrix is applied in the frequency loop, and in the frequency domain, the source has one constant value at each frequency. In Equation (6), if we assume that the waveform of the source is known correctly, the source is finally multiplied by the virtual source ( F v s ) and the residual wave field ( u s - d s ) once. The denominator is multiplied twice each, canceling each other out and removing them efficiently. This means that when source waveform inversion is successfully performed, method 1 eliminates the source's effect so that very low to high frequency data can contribute to the inversion. Therefore, in the previous study, Method 1 provided better results than Method 2, in which the effect of the source was not eliminated.

하지만 방법1의 가장 큰 단점은 앞서 언급했듯이 모든 주파수 성분에 대한 역산 이론에는 맞지 않는다는 것이다. 만약 우리가 방법2에서 송신원의 파형만을 정확하게 제거했다고 가정하면, 식 (7)와 같이 각 주파수 별로 계산된 최대 급경사 방향이 모두 동일한 슈도헤시안 행렬(모든 주파수에 대해서 더해진)에 의해 스케일링 되기 때문에, 각 주파수 성분 별 최대 급경사 방향의 기여도는 다음 식 (8)에 의해 결정된다.The biggest drawback of Method 1, however, is that it does not fit the inversion theory for all frequency components, as mentioned earlier. If we assume that only the source waveform is correctly removed in method 2, since the maximum steep slopes calculated for each frequency are all scaled by the same pseudohessian matrix (added for all frequencies), as shown in equation (7), The contribution of the maximum steep slope to each frequency component is determined by the following equation (8).

Figure 112012015369898-pat00008
(8)
Figure 112012015369898-pat00008
(8)

하지만 방법1에서는 슈도헤시안 행렬이 주파수루프 안에서 적용되기 때문에 각 주파수 별 최대 급경사 방향의 상대적 크기를 결정하는 데에 헤시안 행렬이 영향을 주게 된다. However, in Method 1, since the pseudohessian matrix is applied in the frequency loop, the Hessian matrix influences the relative magnitude of the maximum steepest slope for each frequency.

Figure 112012015369898-pat00009
(9)
Figure 112012015369898-pat00009
(9)

즉, 방법1에서는 식 (9)에 의해 각 주파수 성분 별 최대 급경사 방향의 기여도가 결정되게 되며, 이러한 헤시안 행렬에 의한 가중효과는 역산이론과 반대의 경향을 가진다(식 (8)에서는 가상송신원의 크기에 비례했지만, 식 (9)에서는 가상송신원의 크기에 반비례하는 경향).That is, in the method 1, the contribution of the maximum steep slope for each frequency component is determined by Equation (9), and the weighting effect by the Hessian matrix tends to be opposite to the inversion theory (Equation (8)). (9) tends to be inversely proportional to the size of the virtual source.

이처럼 역산 이론에는 없는 헤시안 행렬에 의한 가중 효과를 제거하기 위해 각 주파수 별로 구해진 최대 급경사 방향 벡터를 최대값으로 정규화하여 사용하는 방법(식 (10))이 더 안정적인 결과를 주고 있다. 이는 공개특허공보 제2010-135602에 개시된 바 있다. In order to remove the weighting effect of the Hessian matrix, which is not included in the inversion theory, the method of normalizing the maximum steep direction vector obtained for each frequency to the maximum value (Eq. (10)) gives more stable results. This has been disclosed in Publication 2010-135602.

Figure 112012015369898-pat00010
(10)
Figure 112012015369898-pat00010
(10)

하지만 각 주파수 별 최대 급경사 방향 벡터를 정규화하여 더하는 것은 다음과 같은 부작용이 있다. 파형 역산 수식에서 중요한 요소 중 하나가 바로 (u s -d s )로 정의되는 잔여파동장이며, 송신원이 제거된 상황을 가정했기 때문에 정확히 말하면, 송신원이 제거된 잔여파동장을 의미한다. 잔여파동장은 실제 지질모델과 현재 이터레이션에서 추정된 모델이 발생시키는 데이터 사이의 차이에 의해 결정되는 값으로, 최대 급경사 방향을 결정하는 수식(식 (5))에서 실제 지질 모델과 역산 모델 사이의 차이를 반영하는 유일한 값이다. 방법1(정규화 방법)은 각 주파수 별 잔여파동장의 세기 또한 평준화시키는 효과가 있으며 이로 인해 실제 지하모델과 추정 중인 모델 사이의 모델링 반응의 차이가 적절하게 반영되지 않는 문제점이 발생한다
However, normalizing the maximum steepness direction vector for each frequency has the following side effects. One of the important factors in the waveform inversion equation is the residual wave field, which is defined as ( u s - d s ), and, precisely, it is the residual wave field from which the source is removed because it is assumed that the source is removed. The residual wave field is a value determined by the difference between the actual geological model and the data generated by the model estimated in the current iteration, and the equation between the actual geological model and the inversion model in equation (5) that determines the maximum steep slope direction. It is the only value that reflects the difference. Method 1 (normalization method) also has the effect of equalizing the residual wave strength for each frequency, which causes a problem that the difference in modeling response between the actual underground model and the estimated model is not properly reflected.

(3)(3) 최대 급경사 방향 분석 : 자코비안 행렬과 잔여 파동장Maximum steep slope analysis: Jacobian matrix and residual wave field

송신원이 제거된 방법2에서, 각 주파수 별 최대 급경사 방향은 다음 식 (11)과 같다. (헤시안 행렬은 모두 동일하게 적용되므로 헤시안 행렬의 주파수 별 가중효과는 고려하지 않는다).In the method 2 in which the transmission source is removed, the maximum steepest inclination direction for each frequency is as shown in Equation (11). (The Hessian matrix is applied equally, so the weighting effect of each frequency of the Hessian matrix is not considered.)

Figure 112012015369898-pat00011
(11)
Figure 112012015369898-pat00011
(11)

위의 식에서 최대 급경사 방향의 상대적 크기는 인위적으로 사용한 송신원이 제거된 자코비안 행렬(J k )과 송신원이 제거된 잔여파동장(u s -d s )에 의해 결정된다. 자코비안 행렬은 해당 주파수와 관련된 단일 주파수의 파동이 전파하는 것과 같은 값의 분포를 가지는데 특징적으로 저주파의 자코비안 행렬은 긴 파장을 가지며, 두꺼운 구조를 감지하기에 유리하며, 고주파의 자코비안 행렬은 짧은 파장을 가지며, 얇은 구조를 감지하기에 유리하다. 따라서 각 주파수 별로 계산된 최대 급경사 방향 또한 그 주파수와 관련된 공간적 분해능을 가지고 있으며, 역산 과정에서 어떤 주파수 성분의 최대 급경사 방향이 더 부각되어야 하는가는 역산 이론에 따라 송신원이 제거된 잔여파동장의 주파수 별 세기의 차이를 반영하여 결정되어야 한다. 그리고 잔여파동장의 주파수 별 세기의 차이는 실제지질모델과 역산 초기모델의 차이에 의해 결정된다.In the above equation, the relative magnitude of the maximum steep slope is determined by the Jacobian matrix J k from which the artificially-used source is removed and the residual wave field u s - d s from which the source is removed. The Jacobian matrix has the same distribution as the propagation of a single frequency wave associated with that frequency. Characteristically, the low-frequency Jacobian matrix has a long wavelength, which is advantageous for detecting thick structures, and the high-frequency Jacobian matrix. Has a short wavelength and is advantageous for detecting thin structures. Therefore, the maximum steep slope calculated for each frequency also has the spatial resolution related to the frequency, and it is the strength of each frequency of the residual wave field from which the source is removed according to the inversion theory. The decision should be made to reflect the differences in. The difference of the intensity of each residual wave field by frequency is determined by the difference between the real geological model and the inversion initial model.

이러한 사실을 확인하기 위해 도 2에 도시한 바와 같은 단순한 모델에 대하여 테스트를 수행하였다.To confirm this fact, a test was performed on a simple model as shown in FIG.

도 2는 본 발명의 지하 매질 구조 추정방법을 테스트하기 위한 단순화된 모델을 도시한 도면이다.2 illustrates a simplified model for testing the method for estimating the underground medium structure of the present invention.

도 2의 (a)는 암염과 같은 두꺼운 지질모델을 단순화시킨 두꺼운 직사각형 모델이고, (b)는 두꺼운 층과 고속도의 얇은 층이 교호하는 지질모델을 단순화시킨 얇은 3층 구조 모델이다. 그리고 역산을 위한 초기 모델은 두 모델의 배경속도와 동일하게 하여, 파동장의 차이가 순수하게 이상체에 의해서만 유발되도록 하였다. FIG. 2 (a) is a thick rectangular model that simplifies a thick geological model such as rock salt, and (b) is a thin three-layer structural model that simplifies a geological model in which a thick layer and a thin layer of high speed are alternated. And the initial model for inversion was made equal to the background velocity of the two models so that the difference in wave field was purely caused by the ideal body.

도 3은 시간 영역에서 송신원이 포함된 잔여파동장을 구한 도면이다.3 is a diagram of a residual wave field including a transmission source in a time domain.

도 3의 (a)는 두꺼운 직사각형 모델에 대해, (b)는 얇은 3층 구조 모델에 대해 얻은 것이다. (a)에서는 두꺼운 이상체에 의해 파장이 긴 저주파의 파가 지배적인 것을 볼 수 있으며, (b)에서는 상대적으로 고주파의 파가 눈에 띄게 부각된다. 그리고 잔여파동장의 세기도 (a)의 경우에 더 강하다. 이러한 차이는 주파수 영역에서도 드러난다.Figure 3 (a) is obtained for the thick rectangular model, (b) is obtained for the thin three-layer structure model. In (a), it can be seen that a low frequency wave with a long wavelength is dominated by a thick ideal body, and in (b), a relatively high frequency wave is noticeable. And the intensity of the residual wave field is stronger in the case of (a). This difference is also apparent in the frequency domain.

도 4는 각 주파수별로 송신원이 제거된 잔여파동장의 제곱으로 정의한 오차를 출력한 도면이다.4 is a diagram illustrating an error defined by the square of the residual wave field from which a transmitter is removed for each frequency.

송신원의 주파수 별 세기 차이를 고려하지 않기 위해 송신원의 영향을 제거하였다. 주파수 별 송신원이 제거된 RMS 오차의 분포는 두꺼운 직사각형 모델의 경우에 큰 차이가 나타난다. 상대적으로 저주파 성분의 오차가 매우 크며, 이를 통해 두꺼운 이상체는 저주파 성분의 오차를 크게 증가시킨다는 것을 확인할 수 있고, 이는 역산 이론에 따라 최대 급경사 방향을 구할 때 저주파 성분의 최대 급경사 방향이 더 가중되어 들어가야 한다는 것을 의미한다. 이와 달리, 얇은 3층 모델의 경우에는 고주파 성분의 오차가 상대적으로 조금 더 크기는 하지만 전체적으로 값이 매우 작고 상대적으로 변화가 적다고 볼 수 있다. 이런 경우, 추가적인 가중 효과가 요구되지 않으며 방법1에서 정규화에 의해 잔여파동장이 평준화되는 것은 역산에 큰 영향을 미치지 않을 것이다.
In order not to consider the difference in intensity by frequency of the transmitter, the influence of the transmitter is removed. The distribution of RMS error with the source removed per frequency is very different for the thick rectangular model. The error of the low frequency component is relatively large, and it can be seen that the thick ideal body greatly increases the error of the low frequency component, which is more weighted when the maximum steep direction of the low frequency component is obtained according to the inversion theory. It means you have to go in. In contrast, in the case of the thin three-layer model, the error of the high frequency component is a little larger, but the overall value is very small and relatively small. In this case, no additional weighting effect is required and leveling the residual wave field by normalization in Method 1 will not have a significant effect on the inversion.

(4) 가중치 기법(4) weighting technique

전술한 바와 같이 송신원이 제거된 잔여파동장의 세기는 지질 모델간의 차이를 반영하여 분포되며, 역산 이론에 따라 각 주파수 별 최대 급경사 방향이 주파수 별 송신원이 제거된 잔여파동장의 세기를 반영해야 함을 확인하였다. 이뿐만 아니라 자코비안 행렬도 주파수에 따라 달라지기 때문에 자코비안 행렬에 의한 가중효과도 고려해주어야 한다. 따라서 본 발명에서는 식 (12)와 같이 정의된 가중 함수를 이용하여 정규화된 각 주파수 별 최대 급경사 방향에 추가적으로 가중치를 주는 방법을 제안한다.As described above, the intensity of the residual wave field from which the source is removed is distributed to reflect the difference between the geological models, and according to the inversion theory, it is confirmed that the maximum steep slope direction for each frequency should reflect the intensity of the residual wave field from which the source is removed for each frequency. It was. In addition, the Jacobian matrix also depends on the frequency, so the weighting effect of the Jacobian matrix should be considered. Therefore, the present invention proposes a method of additionally weighting the maximum steep slope for each frequency normalized using a weighting function defined as Equation (12).

Figure 112012015369898-pat00012
(12)
Figure 112012015369898-pat00012
(12)

vi는 i번째 지점에서의 역전파된 파동장의 복소수 값, 절대값은 진폭, np는 모델링 격자점의 수를 의미한다. 가중치 기법은 기본적으로 송신원 효과를 제거할 수 있는 방법1을 기반으로 하고 있지만, 지하 모델과 역산 과정에서 갱신되고 있는 지질모델간의 차이를 반영하는 가중효과를 주기 위하여 각 최대 급경사 방향 별로 가중치를 주는 방법이다. 이 때 가중치(

Figure 112012015369898-pat00013
)는 위의 식과 같이 송신원이 제거된 역전파 된 파동장 진폭의 평균으로 정의된다. 역전파 된 파동장은 식 (11)에서 마지막 두 개의 항((S -1) T (u s -d s ) )을 의미하며, 물리적으로는 송신원이 제거된 잔여파동장을 송신원으로 하여 역전파 된 파동장을 의미한다. 따라서 역전파 된 파동장의 세기는 송신원이 제거된 잔여파동장의 세기에 비례하며, 가상송신원의 주파수 별 크기 차이를 고려하지 않는다면, 충분히 역산 이론에 근거한 방법2의 가중 효과를 줄 수 있다. v i is the complex value of the back propagated wave field at the i-th point, absolute value is amplitude, and np is the number of modeling grid points. The weighting technique is basically based on Method 1, which can remove the source effect, but the weighting method is applied to each steep slope to give a weighting effect that reflects the difference between the underground model and the geological model being updated during the inversion process. to be. Where the weights (
Figure 112012015369898-pat00013
) Is defined as the average of the backpropagated wave field amplitudes from which the source is removed, as shown above. The back propagated wave field means the last two terms (( S -1 ) T ( u s - d s )) in Equation (11), and is physically back propagated using the residual wave field from which the source is removed. It means wave field. Therefore, the strength of the back propagated wave field is proportional to the strength of the residual wave field from which the transmitter is removed. If the size difference of each frequency of the virtual transmitter is not considered, the weighting effect of the method 2 based on the inversion theory can be sufficiently applied.

가중치 기법의 효율성을 확인하기 위하여 두 가지 방법(방법1, 방법2)과 가중치 기법을 이용하여 도 2에 도시한 두 가지 단순화된 모델에 대해 역산을 수행하였다. In order to verify the efficiency of the weighting technique, two inversions (methods 1 and 2) and weighting techniques were used to perform the inversion on the two simplified models shown in FIG.

도 5는 본 발명의 일실시예에 따른 지하매질 구조 추정 방법의 흐름도이다.5 is a flowchart of a method for estimating an underground medium structure according to an embodiment of the present invention.

도 5를 참조하여 본 발명의 지하매질 구조 추정방법의 흐름을 다시 한번 정리하여 설명하기로 한다.Referring to Figure 5 will be described once again summarized the flow of the method for estimating the structure of the underground medium of the present invention.

다시 말하면, 지하 매질을 통과하여 지표면에서 기록된 탄성파 데이터와 지하 매질에 대한 모델링 데이터를 가지고 파형역산 알고리즘을 적용하여 지하 매질의 구조를 추정하기 위하여, 먼저 측정된 시간 영역의 데이터를 소정의 변환 영역의 데이터로 변환한다(510). 이때 변환 영역은 주파수 도메인, 라플라스 도메인, 또는 라플라스-푸리에 도메인 중 어느 하나가 될 수 있으나 편의상 전술한 예에서는 주파수 도메인으로의 변환을 예를 들어 설명하였다.In other words, in order to estimate the structure of the underground medium by applying the waveform inversion algorithm with the seismic data recorded on the ground surface and the modeling data of the underground medium through the underground medium, the measured time domain data is first converted into a predetermined transform domain. The data is converted into data (510). In this case, the conversion region may be any one of a frequency domain, a Laplace domain, or a Laplace-Fourier domain, but for convenience, the above-described example has been described in terms of the conversion to the frequency domain.

그리고 이렇게 변환된 영역에서의 지하 매질 구조에 대한 측정 데이터와 모델링 데이터의 차이를 이용하여 각 변환변수에서의 목적 함수(objective function)를 정의하고 잔차를 계산한다(520).In operation 520, an objective function of each transformation variable is defined and a residual is calculated using the difference between the measured data and the modeling data of the underground medium structure in the transformed region.

다음으로 변환변수에서의 목적함수를 최소화시키기 위한 최대 급경사 방향을 합산하여 전체 변환변수의 범위에 대한 최대 급경사 방향을 계산하되, 제1가중치를 적용하여 각 변환변수에서의 최대 급경사 방향을 정규화(normalization)하고, 각 변환변수별 최대 급경사 방향에 제2가중치를 추가적으로 적용하여 합산한다(530).Next, the maximum steep slope is calculated by summing the maximum steep directions to minimize the objective function of the transform variables, and normalizing the maximum steep slope directions for each transform variable by applying the first weight. In operation 530, the second weighted value is additionally applied to the maximum steep slope for each conversion variable.

즉, 변환 영역은 주파수 도메인이고, 제1가중치는 상기 지하 매질 구조를 나타내는 소정의 파라미터(pω)에 대한 상기 각 변환변수에서의 최대 급경사 방향 벡터에서 그 절대값이 가장 큰 값의 역수로 정의될 수 있다.That is, the transform region is the frequency domain, and the first weight is defined as the inverse of the largest absolute value of the maximum steepness direction vector in each transform variable for the predetermined parameter p ω representing the underground medium structure. Can be.

또한, 제2가중치는 탄성파 송신원이 제거된 잔여파동장을 송신원으로 하여 역전파된 파동장의 진폭의 평균으로 정의될 수 있다. 전술한 바 있으나, 이를 수식으로 나타내면, 제2가중치(

Figure 112012015369898-pat00014
)는 다음 수학식으로 정의될 수 있다.In addition, the second weighting value may be defined as an average of amplitudes of the back wave propagated wave field using the remaining wave field from which the elastic wave source is removed. As described above, when expressed as a formula, the second weight value (
Figure 112012015369898-pat00014
) Can be defined by the following equation.

Figure 112012015369898-pat00015
Figure 112012015369898-pat00015

여기서 np는 모델링 격자점의 수,

Figure 112012015369898-pat00016
는 역전파된 파동장의 진폭을 나타낸다.Where np is the number of modeling grid points,
Figure 112012015369898-pat00016
Denotes the amplitude of the back propagated wave field.

도 6은 본 발명의 일실시예에 따른 지하매질 구조 추정 장치의 구성도이다.Figure 6 is a block diagram of a structure for estimating the underground medium according to an embodiment of the present invention.

지하 매질을 통과하여 측정된 탄성파 데이터와 상기 지하 매질에 대한 모델링 데이터를 가지고 파형역산 알고리즘을 적용하여 상기 지하 매질의 구조를 추정하는 장치는, 변환부(610), 모델링 데이터 생성부(620), 잔차 산출부(630), 최대 급경사 방향 산출부(640), 갱신부(650)를 포함한다.The apparatus for estimating the structure of the underground medium by applying a waveform inversion algorithm with the seismic data measured through the underground medium and the modeling data for the underground medium includes: a converter 610, a modeling data generator 620, The residual calculator 630 includes a maximum steep slope calculator 640 and an updater 650.

변환부(610)는 측정된 시간 영역의 데이터를 소정의 변환 영역의 데이터로 변환한다.The converter 610 converts the measured data in the time domain into data in the predetermined conversion domain.

모델링 데이터 생성부(620)는 추정하고자 하는 지하 매질을 모델링하기 위해 초기 모델을 설정하고, 변환 영역에서의 모델링 데이터를 산출한다.The modeling data generator 620 sets an initial model to model the underground medium to be estimated, and calculates modeling data in the transform region.

잔차 산출부(630)는 각 변환변수에서의 목적 함수를 정의하고 측정 데이터와 모델링 데이터 사이의 잔차를 산출한다.The residual calculator 630 defines an objective function in each transformation variable and calculates a residual between the measured data and the modeling data.

최대 급경사 방향 산출부(640)은 각 변환변수에서 독립적으로 계산된 최대 급경사 방향을 합산하여 전체 변환변수의 범위에 대한 최대 급경사 방향을 계산하되, 제1가중치를 적용하여 각 변환변수에서의 최대 급경사 방향을 정규화(normalization)하고, 정규화된 각 변환변수별 최대 급경사 방향에 제2가중치를 추가적으로 적용하여 합산하여 최대 급경사 방향을 산출한다.The maximum steep slope calculation unit 640 calculates the maximum steep slope direction for the entire range of the transformation variables by summing the maximum steep slope directions independently calculated from each transformation variable, but applies the first weight to maximize the steep slope slope for each transformation variable. The direction is normalized and the maximum steep slope is calculated by additionally applying a second weighting value to the maximum steep slope for each normalized transform variable.

갱신부(650)는 전체 변환 변수 범위에서 계산된 최대 급경사 방향을 이용하여 목적 함수가 최소가 되는 방향으로 추정하고자 하는 모델 파라미터를 업데이트한다.The updater 650 updates the model parameter to be estimated in the direction in which the objective function is minimum by using the maximum steepness direction calculated in the entire conversion variable range.

도 7은 방법1을 이용해 각각의 모델에 대한 탄성 파형역산을 수행한 결과, 라메상수(λ)에 대해 첫 번째 이터레이션에서 계산된 최대 급경사 방향을 도시한 도면이다.FIG. 7 is a diagram showing the maximum steep slope calculated in the first iteration for the lame constant λ as a result of performing the elastic waveform inversion for each model using Method 1. FIG.

서로 다른 모델임에도 불구하고 최대 급경사 방향 이미지에서 지배적인 주파수 성분이 동일하다. 또한 앞에서 언급했듯이 모델 (b)는 주파수별 송신원이 제거된 RMS 오차의 주파수별 분포가 비교적 평준화되어 있었기 때문에, 정규화에 의한 큰 부작용 없이 비교적 세 개의 얇은 층이 잘 감지가 된다. 하지만 모델 (a) 같은 두꺼운 모델은 저주파 성분의 최대 급경사 방향이 고주파 성분보다 더 가중을 받아야 함에도 불구하고, 얇은 3층 구조모델에서와 동일한 주파수 성분의 최대 급경사 방향이 지배적으로 영향을 주고 있다.Despite the different models, the dominant frequency components in the maximum steepest direction image are the same. Also, as mentioned earlier, in model (b), since the frequency distribution of the RMS error from which the source of frequency was removed is relatively leveled, relatively three thin layers are well detected without significant side effects of normalization. However, in the case of thick models such as model (a), although the maximum steep slope of the low frequency component has to be weighted more than the high frequency component, the steep steep direction of the same frequency component as in the thin three-layer structure model dominates.

도 8은 방법2를 이용해 각각의 모델에 대한 탄성 파형역산을 수행한 결과, 라메상수(λ)에 대해 첫 번째 이터레이션에서 계산된 최대급경사 방향을 도시한 도면이다.FIG. 8 is a diagram showing the maximum steep slope calculated in the first iteration for the lame constant λ as a result of performing the elastic waveform inversion of each model using the method 2. FIG.

방법1의 결과와 마찬가지로 서로 다른 모델임에도 불구하고 최대 급경사 방향 이미지에서 지배적인 주파수가 동일하며 이 주파수는 송신원의 최대주파수의 1/4 정도에 해당한다. 항상 최대주파수의 1/4에 해당하는 최대 급경사 방향이 가중되기 때문에 (b)와 같은 얇은 층 구조를 감지해내지 못하고 있다.As with the results of Method 1, despite the different models, the dominant frequency in the maximum steepness direction image is the same and this frequency corresponds to about 1/4 of the maximum frequency of the source. Since the maximum steep slope corresponding to 1/4 of the maximum frequency is always weighted, a thin layer structure such as (b) cannot be detected.

도 9는 본 발명의 가중치기법을 이용해 각각의 모델에 대한 탄성 파형역산을 수행한 결과, 라메상수(λ)에 대해 첫 번째 이터레이션에서 계산된 최대 급경사 방향을 도시한 도면이다.FIG. 9 is a diagram showing the maximum steep slope calculated in the first iteration for the lame constant λ as a result of performing the elastic waveform inversion of each model using the weighting method of the present invention.

앞의 실제 지질모델과 비교해서 이상체의 특징을 더 잘 잡아내고 있으며, 최대 급경사 방향 이미지에서 지배적인 최대 급경사 방향의 주파수가 상황에 맞게 합리적으로 달라지고 있음을 확인할 수 있다.Compared with the previous geological model, the characteristics of the ideal body are better captured, and it can be seen that the frequency of the steepest steepest direction in the maximum steepness direction image is changed reasonably according to the situation.

도 10은 방법1을 이용해 파형역산을 수행하여 추정한 지하매질에서의 P파 속도를 도시한 도면이다.FIG. 10 is a diagram showing P wave velocity in an underground medium estimated by performing waveform inversion using Method 1. FIG.

앞에서 (b) 모델에 대한 최대 급경사 방향 이미지가 좋았기 때문에 (b) 모델에 대해서는 비교적 성공적으로 역산이 수행되었다. 하지만 (a)와 같은 두꺼운 모델에 대해서, 너무 고주파수 성분의 최대 급경사 방향으로 업데이트 되었기 때문에 그 결과도 이상체의 상부만 물성값을 제대로 찾았다. 그 결과 실제 모델에 비해 더 얇게 역산이 되었고, 이상체의 하부 속도도 왜곡이 심한 것을 확인할 수 있다.The inversion was performed relatively successfully for the model (b) because (b) the image of the maximum steep slope was good for the model. However, for the thick model like (a), because it was updated in the direction of maximum steep slope of high frequency components, only the upper part of the ideal body correctly found the property value. As a result, the inversion was thinner than the actual model, and the lower speed of the ideal body was also severely distorted.

도 11은 방법2를 이용해 파형역산을 수행하여 추정한 지하매질에서의 P파 속도를 도시한 도면이다.FIG. 11 is a diagram showing P wave velocity in an underground medium estimated by performing waveform inversion using Method 2. FIG.

도 11을 참조하면 이 경우는 송신원의 주파수 별 진폭차이에 의한 가중 효과 때문에 두 경우 모두에서 파형역산에 실패했음을 알 수 있다. Referring to FIG. 11, it can be seen that in this case, waveform inversion failed in both cases due to the weighting effect of the amplitude difference for each frequency of the transmission source.

도 12는 가중치기법을 이용해 파형역산을 수행하여 추정한 지하매질에서의 P 속도를 도시한 도면이다.FIG. 12 is a diagram illustrating P velocity in an underground medium estimated by performing waveform inversion using a weighting technique.

앞의 두 가지 방법과 비교해 역산 결과가 훨씬 개선되었음을 확인할 수 있고 고주파 성분의 잡음(artifact)으로 인한 지질구조의 왜곡도 덜한 것을 확인할 수 있다. 이는 앞서 확인했듯이 가중치 기법에서 각 주파수별 최대 급경사 방향의 가중치가 지질모델간의 차이를 반영하여 상대적으로 잘 결정되어 더 개선된 최대 급경사 방향이 구해졌기 때문이다.Compared to the previous two methods, the inversion result is much improved, and the distortion of the geological structure due to the high frequency artifacts is less. This is because, as previously confirmed, in the weighting technique, the weight of the maximum steepness of each frequency is relatively well determined by reflecting the difference between the geological models, and thus the improved maximum steepness of the steepness is obtained.

도 13은 가정한 모델링 영역의 정 중앙(3km)에서 획득한 P파 속도 단면을 도시한 도면이다.FIG. 13 is a diagram showing a P-wave velocity cross section obtained at the center (3 km) of the assumed modeling region.

도 13을 참조하면, 본 발명의 가중치 기법(weighting method)을 이용해 추정한 지하 매질의 탄성파 속도가 실제 모델에 더 잘 부합하는 것을 확인할 수 있다.Referring to FIG. 13, it can be seen that the seismic velocity of the underground medium estimated using the weighting method of the present invention is better suited to the actual model.

도 14는 파형역산이 끝난 후, 최종적으로 송신원이 제거된 RMS 오차를 주파수 별로 출력한 도면이다.FIG. 14 is a view showing, by frequency, an RMS error from which a transmission source is finally removed after waveform inversion is completed. FIG.

도 14를 참조하면, (b) 모델의 경우 고주파 성분의 오차가 진동하는 양상이 심하기 때문에 정확한 비교는 힘들지만, 가중치 기법을 이용하여 파형역산을 수행했을 때 비교적 오차가 더 잘 떨어졌음을 확인할 수 있다. 특히 (a) 모델의 경우에, 기존 방법에서는 두꺼운 이상체에 의해 유발되는 저주파 성분의 오차가 충분히 떨어지지 못한 반면에 가중치 기법에서는 모든 주파수 성분에 대해서 오차가 더 낮은 수준까지 떨어졌고 특히 저주파 성분의 오차가 상대적으로 크게 감소하였다. 이는 가중치 기법으로 구한 최대 급경사 방향이 더 전체 최소값(global minimum)에 가깝다는 것을 의미한다.Referring to FIG. 14, in the case of the model (b), since an error of a high frequency component vibrates severely, accurate comparison is difficult, but when the waveform inversion is performed using a weighting technique, the error is relatively better. In particular, in the case of the model (a), the error of the low frequency component caused by the thick ideal body is not sufficiently reduced in the conventional method, while the weighting technique has fallen to the lower level for all the frequency components, especially the low frequency component. Decreased significantly. This means that the maximum steep slope obtained by the weighting technique is closer to the global minimum.

도 15는 SEG/EAGE 암염모델의 P파 속도와 S파 속도 모델로 탄성 파형역산에서 실제 지질모델로 가정한 모델을 출력한 도면이다.FIG. 15 is a view showing a model assumed as a real geological model in elastic waveform inversion as a P-wave velocity and an S-wave velocity model of the SEG / EAGE rock salt model.

SEG/EAGE 암염모델은 기존의 연구에서 암염구조를 반영하지 못한 부정확한 초기모델을 사용할 경우, 암염의 정확한 속도 값을 찾지 못하고 암염이 더 얇게 역산되고 암염구조 하부의 속도도 부정확하게 추정되는 등 탄성 파형역산이 성공적으로 수행된 사례가 없다. SEG / EAGE rock salt model has elasticity such that when the incorrect initial model that does not reflect the rock salt structure in previous studies does not find the exact rate of rock salt, the rock salt is inverted thinner and the speed of the rock salt structure is incorrectly estimated. No waveform inversion has been successfully performed.

도 16은 초기모델의 P파 속도가 1.5 ­ 3.06 km로 선형적으로 증가하도록 설정했을 때 방법1을 이용해 암염모델에 대한 탄성 파형역산을 수행한 결과 추정한 P파 속도와 S파 속도를 도시한 도면이다.FIG. 16 is a diagram showing P-wave velocity and S-wave velocity estimated by performing elastic waveform inversion on rock salt model using Method 1 when P-wave velocity of initial model is set to increase linearly to 1.5 3.06 km. to be.

방법1을 사용하여 암염모델에 대해 탄성 파형역산을 수행한 결과 기존 연구 사례에서 보고된 바와 같이 암염의 두께가 얇게 역산되고, 하부의 속도 구조도 명확하게 드러나지 않는 것을 확인할 수 있다.As a result of performing the elastic waveform inversion on the rock salt model using the method 1, the thickness of the rock salt was inverted thinly and the velocity structure of the lower part was not clearly revealed, as reported in the previous case.

도 17은 방법1을 이용해 암염모델에 대한 탄성 파형역산을 수행한 결과 파라미터 라메상수(λ와 μ)에 대해 300번째 이터레이션에서 계산된 최대 급경사 방향을 도시한 도면이다.FIG. 17 is a diagram showing the maximum steep slope calculated at the 300 th iteration for the parameter lame constants (λ and μ) as a result of performing elastic waveform inversion on the rock salt model using Method 1. FIG.

도 17을 참조하면 방법1에 의해 계산된 최대 급경사 방향이 두꺼운 암염구조의 정확한 윤곽을 잡아내지 못하고 있으며 이로 인해 역산 결과에서도 도 16과 같이 속도 구조가 부정확하게 추정되는 것을 확인할 수 있다.Referring to FIG. 17, it can be seen that the maximum steep slope calculated by Method 1 does not accurately capture a thick rock salt structure, and thus, the velocity structure is incorrectly estimated as shown in FIG. 16 even in the inversion result.

도 18은 초기모델의 P파 속도가 1.5 ­ 3.06 km로 선형적으로 증가하도록 설정했을 때 가중치기법을 이용해 암염모델에 대한 탄성 파형역산을 수행한 결과 추정한 P파 속도와 S파 속도를 도시한 도면이다.FIG. 18 is a diagram illustrating P-wave velocity and S-wave velocity estimated by performing elastic waveform inversion of rock salt model using weighting method when P-wave velocity of initial model is set to increase linearly to 1.5 3.06 km to be.

방법1을 이용한 탄성 파형역산과 달리 가중치기법을 이용했을 때 암염모델의 속도구조를 더 정확하게 추정하였음을 확인할 수 있다.Unlike the elastic waveform inversion using the method 1, the weight structure is used to estimate the velocity structure of the rock salt model more accurately.

도 19는 가중치기법을 이용해 암염모델에 대한 탄성 파형역산을 수행한 결과 파라미터 라메상수(λ와 μ)에 대해 300번째 이터레이션에서 계산된 최대 급경사 방향을 도시한 도면이다.FIG. 19 is a diagram showing the maximum steep slope calculated at the 300th iteration for the parameter lame constants (λ and μ) as a result of performing the elastic waveform inversion on the rock salt model using the weighting technique.

도 19를 참조하면 가중치기법을 이용하여 계산한 최대 급경사 방향이 실제 암염모델을 더 정확하게 반영하고 있음을 확인할 수 있다.Referring to FIG. 19, it can be seen that the maximum steep slope calculated using the weighting method more accurately reflects the actual rock salt model.

도 20은 각각의 방법을 이용해 암염모델에 대한 탄성 파형역산을 수행한 결과 추정한 6 km와 9 km 지점에서의 P파 속도 단면을 도시한 도면이다.FIG. 20 is a diagram showing P-wave velocity cross-sections estimated at 6 km and 9 km as a result of performing elastic waveform inversion on the rock salt model using the respective methods.

도 20을 참조하면 가중치기법을 이용했을 때 그 동안 파형역산이 어려웠던 암염모델에 대해 지하의 P파 속도 단면을 더 정확하게 추정할 수 있음을 확인할 수 있다.Referring to FIG. 20, it can be seen that the P-wave velocity cross section of the basement can be estimated more accurately for the rock salt model in which waveform inversion was difficult during the weighting technique.

이제까지 본 발명의 바람직한 실시 예들을 중심으로 살펴보았다. 본 발명이 속하는 기술 분야에서 통상의 지식을 가진 자는 본 발명이 본 발명의 본질적인 특성에서 벗어나지 않는 범위에서 변형된 형태로 구현될 수 있음을 이해할 수 있을 것이다. 그러므로 개시된 실시예들은 한정적인 관점이 아니라 설명적인 관점에서 고려되어야 한다. 본 발명의 범위는 전술한 설명이 아니라 특허청구범위에 나타나 있으며, 그와 동등한 범위 내에 있는 모든 차이점은 본 발명에 포함된 것으로 해석되어야 할 것이다.So far I looked at the center of the preferred embodiment of the present invention. It will be understood by those skilled in the art that various changes in form and details may be made therein without departing from the spirit and scope of the invention as defined by the appended claims. Therefore, the disclosed embodiments should be considered in an illustrative rather than a restrictive sense. The scope of the present invention is defined by the appended claims rather than by the foregoing description, and all differences within the scope of equivalents thereof should be construed as being included in the present invention.

Claims (10)

지하 매질을 통과하여 측정된 탄성파 데이터와 상기 지하 매질에 대한 모델링 데이터를 가지고 파형역산 알고리즘을 적용하여 상기 지하 매질의 구조를 추정하는 방법에 있어서,
상기 측정된 탄성파 데이터를 시간 영역에서 소정의 변환 영역의 데이터로 변환하고, 상기 변환 영역에서의 상기 지하 매질 구조에 대한 측정 데이터와 모델링 데이터의 차이를 이용하여 각 변환변수에서의 목적 함수(object function)를 정의하고 잔차를 계산하는 단계; 및
상기 각 변환변수에서의 목적 함수를 최소화하기 위해 계산된 최대 급경사 방향을 합산하여 전체 변환변수의 범위에 대한 최대 급경사 방향을 계산하되, 제1가중치를 적용하여 상기 각 변환변수에서의 최대 급경사 방향을 정규화(normalization)하고, 상기 정규화된 각 변환변수별 최대 급경사 방향에 제2가중치를 추가적으로 적용하여 합산하는 단계를 포함하는 것을 특징으로 하는 지하 매질 구조를 추정하는 방법.
In the method for estimating the structure of the underground medium by applying a waveform inversion algorithm with the seismic data measured through the underground medium and the modeling data for the underground medium,
Converts the measured acoustic wave data into data of a predetermined transformation region in a time domain, and uses an object function of each transformation variable by using a difference between measured data and modeling data for the underground medium structure in the transformation region. ) And calculating the residuals; And
In order to minimize the objective function in each of the conversion variables, the calculated maximum steep slope is calculated by summing the maximum steep slope directions for the entire range of the transformation variables. And normalizing and adding the second weighting value to the maximum steep slope for each normalized conversion variable.
제1항에 있어서,
상기 변환 영역은 주파수 도메인, 라플라스 도메인, 또는 라플라스-푸리에 도메인을 포함하는 것을 특징으로 하는 지하 매질 구조를 추정하는 방법.
The method of claim 1,
And wherein said transform region comprises a frequency domain, a Laplace domain, or a Laplace-Fourier domain.
제2항에 있어서,
상기 변환 영역은 주파수 도메인이고, 상기 제1가중치는 상기 지하 매질 구조를 나타내는 소정의 파라미터(pω)에 대한 상기 각 변환변수에서 계산된 최대 급경사 방향 벡터에서 그 절대값이 제일 큰 값의 역수로 정의되는 것을 특징으로 하는 지하 매질 구조를 추정하는 방법.
3. The method of claim 2,
The transform region is in the frequency domain, and the first weight value is an inverse of the absolute value of the largest absolute value in the maximum steepness direction vector calculated from the respective transform variables for the predetermined parameter p ω representing the underground medium structure. Method for estimating the underground medium structure, characterized in that it is defined.
제2항에 있어서,
상기 변환 영역은 주파수 도메인이고, 상기 제2가중치는 상기 탄성파 송신원이 제거된 측정 데이터와 모델링 데이터 사이의 잔여파동장을 송신원으로 하여 계산된 역전파된 파동장 진폭의 평균으로 정의되는 것을 특징으로 하는 지하 매질 구조를 추정하는 방법.
3. The method of claim 2,
The transform region is a frequency domain, and the second weighting value is defined as an average of backpropagated wave field amplitudes calculated using a residual wave field between the measured data from which the elastic wave transmitter is removed and modeling data as a transmission source. How to estimate the structure of the underground medium.
제4항에 있어서,
상기 제2가중치(
Figure 112012015369898-pat00017
)는 다음 수학식으로 정의되는 것을 특징으로 하는 지하 매질 구조를 추정하는 방법.
Figure 112012015369898-pat00018

여기서 np는 모델링 격자점의 수,
Figure 112012015369898-pat00019
는 역전파된 파동장의 진폭을 나타낸다.
5. The method of claim 4,
The second weight value (
Figure 112012015369898-pat00017
) Is a method for estimating the underground medium structure, characterized in that defined by the following equation.
Figure 112012015369898-pat00018

Where np is the number of modeling grid points,
Figure 112012015369898-pat00019
Denotes the amplitude of the back propagated wave field.
지하 매질을 통과하여 지표면에서 기록된 탄성파 데이터와, 가정한 초기 지하 매질에 대한 모델링 데이터를 가지고 파형역산 알고리즘을 적용하여 상기 지하 매질의 구조를 추정하는 장치에 있어서,
상기 측정된 탄성파 데이터를 시간 영역에서 소정의 변환 영역의 데이터로 변환하는 변환부;
상기 추정하고자 하는 지하 매질을 모델링하는 초기 모델을 설정하고, 상기 변환 영역에서의 모델링 데이터를 산출하는 모델링 데이터 생성부;
상기 변환 영역의 각 변환변수에서의 목적함수를 정의하고 상기 측정된 탄성파 데이터와 상기 모델링 데이터간의 잔차를 계산하는 잔차 산출부;
상기 각 변환변수에서의 최대 급경사 방향을 독립적으로 계산한 뒤 합산하여 전체 변환변수의 범위에 대한 최대 급경사 방향을 계산하되, 제1가중치를 적용하여 상기 각 변환변수에서의 최대 급경사 방향을 정규화(normalization)하고, 상기 각 정규화된 각 변환변수별 최대 급경사 방향에 제2가중치를 추가적으로 적용한 뒤 합산하여 최대 급경사 방향을 계산하는 최대 급경사 방향 산출부; 및
상기 전체 변환 변수 범위에서 계산된 최대 급경사 방향을 이용하여 상기 목적 함수가 최소가 되는 방향으로 상기 모델링 데이터를 업데이트하는 갱신부를 포함하는 것을 특징으로 하는 지하 매질 구조를 추정하는 장치.
In the device for estimating the structure of the underground medium by applying a waveform inversion algorithm with the seismic data recorded on the ground surface through the underground medium and the modeling data for the assumed initial underground medium,
A converting unit converting the measured acoustic wave data into data of a predetermined conversion region in a time domain;
A modeling data generator configured to set an initial model for modeling the underground medium to be estimated and to calculate modeling data in the transformed region;
A residual calculation unit defining an objective function in each transformation variable of the transformation region and calculating a residual between the measured acoustic wave data and the modeling data;
Calculate the maximum steep slope for the entire range of the transformed variable by independently calculating and summing up the maximum steep slope in each of the transform variables, and normalizing the maximum steep slope in each of the transform variables by applying a first weight value. A maximum steep direction calculation unit configured to calculate a maximum steep slope direction by adding a second weighting value to the maximum steep slope direction for each normalized conversion variable and adding the sum; And
And an updater configured to update the modeling data in a direction in which the objective function is minimized by using the maximum steep slope calculated in the entire transformation variable range.
제6항에 있어서,
상기 변환 영역은 주파수 도메인, 라플라스 도메인, 또는 라플라스-푸리에 도메인을 포함하는 것을 특징으로 하는 지하 매질 구조를 추정하는 장치.
The method according to claim 6,
Wherein the transform region comprises a frequency domain, a Laplace domain, or a Laplace-Fourier domain.
제7항에 있어서,
상기 변환 영역은 주파수 도메인이고, 상기 제1가중치는 상기 지하 매질 구조를 나타내는 소정의 파라미터(pω)에 대한 상기 각 변환변수에서 계산된 최대 급경사 방향 벡터에서 그 절대값이 제일 큰 값의 역수로 정의되는 것을 특징으로 하는 지하 매질 구조를 추정하는 장치.
The method of claim 7, wherein
The transform region is in the frequency domain, and the first weight value is an inverse of the absolute value of the largest absolute value in the maximum steepness direction vector calculated from the respective transform variables for the predetermined parameter p ω representing the underground medium structure. Apparatus for estimating the underground medium structure, characterized in that it is defined.
제7항에 있어서,
상기 변환 영역은 주파수 도메인이고, 상기 제2가중치는 상기 탄성파 송신원이 제거된 측정 데이터와 모델링 데이터 사이의 잔여 파동장을 송신원으로 하여 계산된 역전파된 파동장 진폭의 평균으로 정의되는 것을 특징으로 하는 지하 매질 구조를 추정하는 장치.
The method of claim 7, wherein
The transform region is a frequency domain, and the second weight value is defined as an average of backpropagated wave field amplitudes calculated using a residual wave field between the measured data from which the elastic wave transmitter is removed and modeling data as a transmission source. Device for estimating the structure of the underground medium.
제9항에 있어서,
상기 제2가중치(
Figure 112013043079146-pat00020
)는 다음 수학식으로 정의되는 것을 특징으로 하는 지하 매질 구조를 추정하는 장치.
Figure 112013043079146-pat00021

여기서 np는 모델링 격자점의 수,
Figure 112013043079146-pat00022
는 역전파된 파동장의 진폭을 나타낸다.
10. The method of claim 9,
The second weight value (
Figure 112013043079146-pat00020
) Is a device for estimating the underground medium structure, characterized in that defined by the following equation.
Figure 112013043079146-pat00021

Where np is the number of modeling grid points,
Figure 112013043079146-pat00022
Denotes the amplitude of the back propagated wave field.
KR1020120019214A 2012-02-24 2012-02-24 Method and apparatus of estimating underground structure using a plurality of weighting value KR101318994B1 (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
KR1020120019214A KR101318994B1 (en) 2012-02-24 2012-02-24 Method and apparatus of estimating underground structure using a plurality of weighting value

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
KR1020120019214A KR101318994B1 (en) 2012-02-24 2012-02-24 Method and apparatus of estimating underground structure using a plurality of weighting value

Publications (2)

Publication Number Publication Date
KR20130097504A KR20130097504A (en) 2013-09-03
KR101318994B1 true KR101318994B1 (en) 2013-10-16

Family

ID=49449857

Family Applications (1)

Application Number Title Priority Date Filing Date
KR1020120019214A KR101318994B1 (en) 2012-02-24 2012-02-24 Method and apparatus of estimating underground structure using a plurality of weighting value

Country Status (1)

Country Link
KR (1) KR101318994B1 (en)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9910174B2 (en) 2014-07-25 2018-03-06 Seoul National University R&Db Foundation Seismic imaging apparatus and method for performing iterative application of direct waveform inversion
CN106597484B (en) * 2016-12-12 2018-11-23 武汉大学 The precise quantification method that thermal expansion effects influence GPS coordinate time series

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6999880B2 (en) 2003-03-18 2006-02-14 The Regents Of The University Of California Source-independent full waveform inversion of seismic data
KR20080114488A (en) * 2007-06-26 2008-12-31 신창수 An apparatus for imaging a subsurface structure using waveform inversion in the laplace domain and methods thereof
KR20100135602A (en) * 2009-06-17 2010-12-27 서울대학교산학협력단 Apparatus and method for imaging a subsurface using waveform inversion

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6999880B2 (en) 2003-03-18 2006-02-14 The Regents Of The University Of California Source-independent full waveform inversion of seismic data
KR20080114488A (en) * 2007-06-26 2008-12-31 신창수 An apparatus for imaging a subsurface structure using waveform inversion in the laplace domain and methods thereof
KR20100135602A (en) * 2009-06-17 2010-12-27 서울대학교산학협력단 Apparatus and method for imaging a subsurface using waveform inversion

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
논문.2011 *

Also Published As

Publication number Publication date
KR20130097504A (en) 2013-09-03

Similar Documents

Publication Publication Date Title
KR101092668B1 (en) Apparatus and method for imaging a subsurface using waveform inversion
US9864083B2 (en) Beat tone full waveform inversion
US9829592B2 (en) Seismic imaging with visco-acoustic reverse-time migration using pseudo-analytical method
Chen et al. Detecting a known near-surface target through application of frequency-dependent traveltime tomography and full-waveform inversion to P-and SH-wave seismic refraction data
EP3586169B1 (en) Generating geophysical images using directional oriented wavefield imaging
US11215726B2 (en) Inversion with exponentially encoded seismic data
US20150293245A1 (en) Method and device for the generation and application of anisotropic elastic parameters in horizontal transverse isotropic (hti) media
US20180164453A1 (en) Method for Improved Geophysical Investigation
US20170184748A1 (en) A method and a computing system for seismic imaging a geological formation
KR20120019055A (en) Apparatus and method for imaging a subsurface using frequency domain reverse time migration in an elastic medium
EA032186B1 (en) Seismic adaptive focusing
US11733413B2 (en) Method and system for super resolution least-squares reverse time migration
RU2570827C2 (en) Hybrid method for full-waveform inversion using simultaneous and sequential source method
US11635540B2 (en) Methods and devices performing adaptive quadratic Wasserstein full-waveform inversion
KR101318994B1 (en) Method and apparatus of estimating underground structure using a plurality of weighting value
US9557429B2 (en) Systems and methods for reducing noise in a seismic vibratory source
US8391564B2 (en) Providing an imaging operator for imaging a subterranean structure
Liu et al. Absolute acoustic impedance inversion using convolutional neural networks with transfer learning
US9348048B2 (en) Seismic data processing and apparatus
CN115170428A (en) Noise reduction method for acoustic wave remote detection imaging graph
KR101355107B1 (en) Method and apparatus of estimating underground structure using denoising
US9229121B2 (en) Seismic imaging system for acoustic-elastic coupled media using accumulated Laplace gradient direction
KR101352621B1 (en) seismic imaging method considering a contour of the sea bottom
US20120323541A1 (en) Seismic imaging method considering a contour of the sea bottom
WO2015124960A1 (en) Systems and methods for improved inversion analysis of seismic data

Legal Events

Date Code Title Description
A201 Request for examination
E902 Notification of reason for refusal
E701 Decision to grant or registration of patent right
GRNT Written decision to grant
FPAY Annual fee payment

Payment date: 20160217

Year of fee payment: 4

FPAY Annual fee payment

Payment date: 20170925

Year of fee payment: 5

FPAY Annual fee payment

Payment date: 20181002

Year of fee payment: 6