KR101352621B1

KR101352621B1 - 해저면 지형을 반영한 지하구조 영상화 방법

Info

Publication number: KR101352621B1
Application number: KR1020120063898A
Authority: KR
Inventors: 신창수
Original assignee: 서울대학교산학협력단
Priority date: 2011-06-14
Filing date: 2012-06-14
Publication date: 2014-01-17
Also published as: KR20120138702A

Abstract

본 발명은 송신원으로부터 출력된 파동이 지하 구조를 전파한 후 수신된 측정 데이터를 처리하여 지하 구조를 영상화하는 탄성파 탐사(Seismic imaging) 기술에 관련된다. 일 양상에 따르면, 파동 방정식의 계수 행렬은 해저면 부근의 광역 격자(Grid) 중 적어도 일부 내에서 해저면의 윤곽을 반영하여 산출된다. 격자 내부의 보다 섬세한 해저면 윤곽을 반영하기 때문에 보다 정확하게 해저면에서 반사되거나 투과한 신호들을 정확히 추정할 수 있다. 나아가 계산 부하도 최소화될 수 있다.

Description

해저면 지형을 반영한 지하구조 영상화 방법{seismic imaging method considering a contour of the sea bottom}

본 발명은 송신원으로부터 출력된 파동이 지하 구조를 전파한 후 수신된 측정 데이터를 처리하여 지하 구조를 영상화하는 탄성파 탐사(Seismic imaging) 기술에 관련된다.

파형 역산에 의한 지하구조의 영상화 기술이 알려져 있다. 예를 들어 본 발명자에 의해 발명되어 한국 특허청에 2009. 6. 17.자로 출원되어 2011. 12. 5.자로 특허제1,092,668호로 등록된 발명이 알려져 있다. 이 발명은 또한 미국 특허청에 출원번호 12/817,799호로 출원되어 있다.

이에 따르면 송신원으로부터 송출된 저주파 신호가 지하 구조를 통과한 후 반사된 반사파를 하이드로폰 어레이와 같은 수신기를 통해 측정한 측정 데이터를 사용해 지하 구조의 모델링 파라메터를 구한다. 파동 방정식의 계수들은 이 파동이 전파되는 지하 매질의 밀도와 같은 모델링 파라메터들로 구성된다. 파동 방정식의 모델링 파라메터들이 파형 역산에 의해 계산된다. 파형 역산에 의하면, 파동방정식의 해인 모델링 데이터와 측정 데이터간의 차이에 관한 잔차 함수(residual function)를 최소화하는 방향으로 모델링 파라메터들이 반복적으로(iteratively) 갱신됨에 의해 산출된다.

위 개시된 발명에 따르면, 모델링 파라메터로부터 구해지는 계수 행렬을 가진 파동 방정식의 해인 모델링 데이터와 수신기에서 측정된 측정 데이터와의 차이에 관한 잔차 함수를 최소화하는 방향으로 모델링 파라메터의 갱신을 반복하여 모델링 파라메터가 구해진다. 모델링 파라메터를 구하기 위해 먼저 모델링 파라메터로부터 파동 방정식의 계수 행렬이 계산된다. 이후에 이 계수 행렬을 가진 파동 방정식을 주어진 송신원 정보로 풀어서 그 해인 모델링 데이터를 구한다. 다음으로, 측정 데이터와 모델링 데이터 계산 단계에서 계산된 모델링 데이터와의 잔차에 관한 잔차 함수를 구한다. 만약 산출된 잔차 함수의 값이 기준 이상이면, 현재 파동 방정식의 모델링 파라메터를 잔차 함수를 최소화하는 방향으로 갱신한다. 만약 산출된 잔차 함수의 값이 기준 이하이면 현재 모델링 파라메터를 출력값으로 출력한다.

흔히 이러한 파형 역산은 광역 격자(global grid)라고 불리는 격자점들을 기준으로 수행된다. 종래 해저면은 이 거친 격자점들에 맞추어 모델링되었다. 따라서 해저면에서 반사되거나 그 해저면을 투과한 신호들을 정확히 추정하는 것이 어려웠다.

제시된 발명은 계산량의 부담을 크게 늘이지 않으면서도 해저면의 윤곽을 정확히 반영하여 해저면을 투과한 파동장을 정확히 추정하는 것을 목적으로 한다.

이러한 목적을 달성하기 위한 일 양상에 따르면, 파동 방정식의 계수 행렬은 해저면 부근의 광역 격자(Grid) 중 적어도 일부 내에서 해저면의 윤곽을 반영하여 산출된다. 격자 내부의 보다 섬세한 해저면 윤곽을 반영하기 때문에 보다 정확하게 해저면에서 반사되거나 투과한 신호들을 정확히 추정할 수 있다. 나아가 계산 부하도 최소화될 수 있다.

구체적인 양상에 따르면, 계수 행렬의 산출은 질량 행렬의 산출을 포함하는데, 질량 행렬은 하나의 광역 격자를 해저면 상부의 매질인 제 1 영역과 하부의 매질인 제 2 영역으로 나누고, 두 영역에 대해 수치적분법(numerical integration method)을 적용하여 산출된다.

또다른 양상에 따르면, 수치적분법은 가우스 구적법(Gaussain Quadrature Integration Method)일 수 있다.

본 발명에 따르면, 격자 내부의 보다 섬세한 해저면 윤곽을 반영하기 때문에 보다 정확하게 해저면에서 반사되거나 투과한 신호들을 정확히 추정할 수 있다. 나아가 계산 부하도 최소화될 수 있다.

도 1은 일 실시예에 따른 지하구조 영상화 방법의 개략적인 구성을 도시한 흐름도이다.
도 2는 해저면에 의해 구분된 직사각형 공간요소(cubic elements)의 2차원 단면도를 예시한다.

전술한, 그리고 추가적인 본 발명의 양상들은 후술하는 실시예들을 통해 더욱 명확해질 것이다. 이하에서는 이러한 양상들을 첨부된 도면을 참조하여 기술되는 실시예들을 통해 상세히 설명하기로 한다.

도 1은 일 실시예에 따른 지하구조 영상화 방법의 개략적인 구성을 도시한 흐름도이다. 일 실시예에 따른 지하구조 영상화 방법은 파형 역산 과정을 포함한다. 일 양상에 따르면, 파형 역산은 모델링 파라메터로부터 구해지는 계수 행렬을 가진 파동 방정식의 해인 모델링 데이터와 수신기에서 측정된 측정 데이터와의 차이에 관한 잔차 함수를 최소화하는 방향으로 모델링 파라메터의 갱신을 반복하여 모델링 파라메터를 구한다.

Shin, C. S., & Cha, Y. H., 2008. Waveform inversion in the Laplace domain, Geophys. J. Int., 173, 922-931.에는 라플라스 영역의 파형 역산을 소개하고 있다. 본 발명은 이에 한정되지 않으며, 이 실시예는 예시적인 것이다. 시간 영역에서 라플라스 변환된 파동장은 다음과 같이 표현된다.

(1)

여기서

는 라플라스 영역에서의 파동장, u(t)는 시간 영역에서의 파동장, t 는 시간, s 는 라플라스 감쇄 정수이다.

시간 영역의 파동 방정식을 라플라스 영역으로 변환하여 라플라스 영역에서의 파동 방정식을 구할 수 있다.

(2)

여기서 c(x,y,z) 는 P-파 속도이고, u(x,y,z,t)는 압력장(pressure field)이며, f(x,y,z,t)는 송신원 함수이고, 각 글자 위의 모자(hat) 표시는 라플라스 영역으로 변환된 변수임을 나타낸다.

라플라스 영역의 파동 방정식은 유한요소법(finite element method)에 의해 풀 수 있다. 가중 잔차법(weighted residual method)을 적용하여 파동 방정식에 등가인 수정된 공식을 유도한다. 가중 잔차법을 식(2)에 적용하기 위해 잔차 함수를 다음과 같이 정의한다.

(3)

여기서 ▽는 라플라스 연산자를 의미한다.

식 (3)에 임의의 가중 함수 v 를 곱하여, 유연형(weak form)으로 변경하고, 주어진 영역 Ω 에서 적분하면,

(4)

식 (4)를 부분적분하고 자연 경계 조건(natural boundary condition)을 적용하면,

(5)

라플라스 영역으로 변환된 파동장

와, v 는 갈러킨 근사법(Galerkin approximation)을 이용하여 가중 함수 α_j(s)와 β_j(s), 그리고 기저 함수(basis function)들인 φ_j(x,y,z)와 φ_i(x,y,z)의 가중합으로 다음과 같이 나타낼 수 있다.

(6)

여기서 임의 함수인 v=1로 가정하고, 식(6)을 식(5)에 대입하고 재배열하면,

(7)

로 구할 수 있다.

기저 함수의 계수들인 α_j(s)가 벡터이므로, 기본적으로 이 계수들은 파동장을 나타내므로, 식 (7)을 행렬의 곱의 형태로 다음과 같이 표현할 수 있다.

(8)

여기서,

(9)

위 식(9)에서, M은 질량 행렬(mass matrix)이고 K는 경도 행렬(stiffness matrix)이다. 라플라스 영역에서의 파동장은 식 (8)을 다음의 식 (10)에 기술된 바와 같이 풀어서 구할 수 있다.

(10)

도 1은 일 실시예에 따른 지하구조 영상화 방법의 개략적인 구성을 도시한 흐름도이다. 특허제1,092,668호에 기술된 바와 같이, 모델링 파라메터로부터 구해지는 계수 행렬을 가진 파동 방정식의 해인 모델링 데이터와 수신기에서 측정된 측정 데이터와의 차이에 관한 잔차 함수 r을 최소화하는 방향으로 모델링 파라메터의 갱신을 반복하여 모델링 파라메터가 구해진다. 도 1에 도시된 바와 같이, 모델링 파라메터를 구하기 위해 먼저 모델링 파라메터로부터 파동 방정식의 계수 행렬이 계산된다(단계 100~300). 이후에 이 계수 행렬을 가진 파동 방정식을 주어진 송신원 정보로 풀어서 그 해인 모델링 데이터를 구한다(단계 400). 다음으로, 측정 데이터와 모델링 데이터 계산 단계에서 계산된 모델링 데이터와의 잔차에 관한 잔차 함수를 구한다(단계 500).

잔차 함수의 계산에 대해서는 Pyun, S. J., Shin, C. S. & Bednar, J. B., 2007. Comparison of waveform inversion, part3: amplitude approach, Geophys. Prospect., 55, 465-475.에 설명되어 있다. 로그 잔차함수의 적용에 대해서는 Shin, C. S., & Min, D. J., 2006. Waveform inversion using a logarithmic wavefield: Geophysics, 71, R31-R42.에 설명되어 있다.

이후에 산출된 잔차 함수의 값을 기준 값인 R_ref 와 비교한다(단계 600). 만약 산출된 잔차 함수의 값이 기준값인 R_ref 이상이면, 현재 파동 방정식의 모델링 파라메터를 잔차 함수를 최소화하는 방향으로 갱신한다(단계 700). 잔차 함수를 최소화하는 방향을 결정하기 위해, 잔차 함수의 그래디언트(gradient) 방향이 산출된다. 가우스-뉴톤법(Gauss-Newton method), 마쿼트-레벤버그법(Marquardt-Levenverg method), 최급감법(the steepest decent method) 등의 최소자승법(least-square method) 알고리즘들이 적용될 수 있다. 그래디언트 방향을 효율적으로 계산하기 위해 후방 전파법(back-propagation algorithm)을 사용하는 기술이 전술한 Shin & Min 2006 논문에 설명되어 있다. 갱신된 모델링 파라메터를 사용하여 다시 파동 방정식의 계수 행렬이 계산된다(단계 200). 이와 같은 과정들이 반복되다, 만약 산출된 잔차 함수의 값이 기준값인 R_ref 이하이면 현재 모델링 파라메터를 출력값으로 출력한다(단계 800).

일 양상에 따르면, 계수 행렬은 해저면 부근의 광역 격자(Grid) 중 적어도 일부 내에서 해저면의 윤곽을 반영하여 구하는 질량 행렬로부터 계산된다. 보다 구체적인 양상에 따르면, 질량 행렬은 하나의 광역 격자를 해저면 상부의 매질인 제 1 영역과 하부의 매질인 제 2 영역으로 나누고, 두 영역에 대해 수치적분법(numerical integration method)을 적용하여 산출한다.

도 2는 해저면에 의해 구분된 직사각형 공간요소(cubic elements)의 2차원 단면도를 예시한다. 각 공간 요소들은 광역 격자(global grid)에 의해 구획된다. 경사지게 그려진 선, 즉 3개의 사각형 점들을 연결하는 인터페이스는 가정된 해저면을 표시한다. 이 선들이 확장된 수치 적분점들(extended numerical integration points)을 2개의 상이한 그룹 Ω1 및 Ω2로 나눈다. 공간요소 질량 행렬(element mass matrix)은 Ω1 및 Ω2의 각 그룹에 할당된 상이한 속도 모델을 사용하는 수치 적분법에 의해 계산될 수 있다.

해저 경계면 상부의 공간요소에 있어서 매질들은 균일하기 때문에 이러한 구획된 수치적분법이 적용될 필요가 없다. 즉, 해저 경계면 상부와, 탐사선 하부의 공간에 속하는 공간 요소들은 모두 바다이므로 여기서 모델링 파라메터, 예를 들면 밀도 혹은 파동 전파 속도는 균일한 것으로 가정된다. 따라서 적어도 해저 경계면에 대응되는 공간 요소들 중 일부, 특히 경계면이 경사진 공간 요소들에 있어서 본 발명에 따른 구획된 수치적분법이 적용된다. 이에 의해 매 반복처리(iteration) 때 마다 계산되어야 하는 확장된 수치 적분법이 적용되는 공간 요소들의 수가 크게 줄어들고, 이에 의해 전체 지하구조 영상화의 계산 부하를 별로 늘이지 않으면서도 경계면에 의한 오차를 크게 줄일 수 있다. 특히 3차원 지하구조 영상화의 경우 계산량이 매우 크고 해저면 형상의 영향도 크기 때문에 이러한 점은 더욱 중요하다.

또다른 양상에 따르면, 수치적분법은 가우스 구적법(Gaussain Quadrature Integration Method)일 수 있다. 가우스 구적법은 (2n+1)차 임의 함수의 일차원 적분을 n 개의 적분점 좌표와 그 가중치값의 가중선형합(linear combination)dmfh 표현하는 수치적분법이다.

일 양상에 따르면, 해저면을 따라 존재하는 공간 요소들에게 라플라스 영역 모델링과 파형 역산 알고리즘을 적용함에 있어서 임피던스 행렬을 구성하는 질량 행렬을 계산할 때 가우스 구적법이 적용된다. 가우스 구적법에서, 식 (9)의 질량 행렬은 다음과 같이 표현된다.

(11)

식 (11)에 있어서, M_ij ^e 는 공간 요소 질량 행렬이고, p, q, r은 3차원 매질의 한 공간 요소 내에서 가우스 구적분점(gaussian quadrature point)의 인덱스이며, h 는 격자 간격이다. φ_i, φ_j 는 각각 공간 내의 i, j번째 점에서 형상 함수들(shape functions)의 값이다. 각 형상 함수는 하나의 격자점에서 '1'의 값을 갖고 나머지 격자점에서는 모두 '0'의 값을 갖는 함수이다. 각 형상 함수들이 '1'의 값을 갖는 점은 모두 다르다. 적분점들의 국부 좌표는 ξ_p, η_q, ζ_r이고, F(ξ_p,η_q,ζ_r)는 국부 좌표에서의 형상 함수의 값들의 곱이다. 매질 파라메터인 속도 c 가 공간의 함수이므로, 격자 간격이 넓어서 해저면 경계가 그 공간 격자를 통과하는 경우에는 해저면 부근의 공간 요소 내에서 일정한 값을 가지지 않는다. 그러나 통상적인 3차원 라플라스 영역 모델링 기술에서는 하나의 공간 요소 내에서는 동일한 속도를 가지는 것으로 기술하기 때문에 해상도 문제가 야기된다.

하나의 공간 요소에 두 개의 속도를 반영하기 위해서, 해저 경계면 부근에서 질량 행렬을 적분할 때 다음과 같이 각 가우스 구적분점 (ξ_p,η_q,ζ_r)의 속도에 대응되는 F(ξ_p,η_q,ζ_r) 값을 다르게 사용한다.

(12)

여기서 Ω1, Ω2는 해저면에 의해 나뉘어진 공간 요소를 포함하는 영역이다. 이 방법은 속도 성분의 가중치로 속도의 공간적 분포를 사용하여 가중하는 방식으로 해석될 수 있다. 해저면에서 공간 요소 내에서 단 하나의 성분값만을 적용한다면, 수중 속도를 선택하든 지하 속도를 선택하든, 질량 행렬의 값은 이 두 극단 중의 하나가 된다. 따라서 이 두 속도 값을 반영한 중도 값을 사용하는 것이 바람직하다.

이상에서 본 발명은 첨부된 도면을 참조하여 기술되는 실시예를 중심으로 기재되었으나, 이에 한정되는 것은 아니며, 이로부터 자명하게 도출 가능한 많은 변형예들을 포괄하도록 해석된다. 첨부된 청구범위는 이러한 변형예를 포괄하도록 의도되었다.

Claims

모델링 파라메터로부터 구해지는 계수 행렬을 가진 파동 방정식의 해인 모델링 데이터와 수신기에서 측정된 측정 데이터와의 차이에 관한 잔차 함수를 최소화하는 방향으로 모델링 파라메터의 갱신을 반복하여 모델링 파라메터를 구하는 파형 역산 단계를 포함하는 지하구조 영상화 방법에 있어서,
해저면 부근의 광역 격자(Grid) 중 적어도 일부 내에서 해저면의 윤곽을 반영하여 광역 격자를 해저면을 기준으로 해저면 상부의 매질인 제 1 영역과 그 하부의 매질인 제 2 영역으로 나누고, 제 1 영역 및 제 2 영역에 대해 각각 수치적분법(numerical integration method)을 적용하여 산출하는 질량 행렬로부터 계수 행렬을 구하는 단계와;
파동 방정식을 주어진 송신원 정보로 풀어서 그 해인 모델링 데이터를 구하는 모델링 데이터 계산 단계와;
수진기에서 측정된 측정 데이터와 모델링 데이터 계산 단계에서 계산된 모델링 데이터와의 잔차에 관한 잔차 함수를 구하는 잔차 함수 계산 단계와;
잔차 함수 계산 단계에서 산출된 잔차 함수의 값이 기준 이상이면, 현재 파동 방정식의 모델링 파라메터를 잔차 함수를 최소화하는 방향으로 갱신하고, 이하이면 현재 산출된 질량 행렬로부터 지하 구조의 영상화 데이터의 출력값을 산출하여 출력하는 모델링 파라메터 계산 단계;
를 포함하는 지하구조 영상화 방법.
삭제
제 1 항에 있어서, 수치적분법이 가우스 구적법(Gaussain Quadrature Integration Method)인 지하구조 영상화 방법.
제 1 항에 따른 3차원 지하구조 영상화 방법.
삭제
제 4 항에 있어서, 수치적분법이 가우스 구적법(Gaussain Quadrature Integration Method)인 지하구조 영상화 방법.
삭제
제 1, 3, 6 항 중 어느 한 항에 따른 방법이 구현된 컴퓨터로 판독 가능한 프로그램이 저장된 기록 매체.