아이코날 방정식

1. 개요[편집]

아이코날 방정식 Eikonal Equation
형태	1계 비선형 쌍곡형 PDE (해밀턴-야코비 계열)
미지수	도달 시간 T(x) — 파면이 언제 오는가
출처	파동방정식의 고주파(WKB) 극한
해의 개념	점성해 (크랜들·라이온스, 1983)
대표 해법	고속행진법 O(N log N) · 고속 스위핑법 O(N)
주 응용	초동주시 계산 · 경로계획 · 부호거리함수

파동방정식을 다 풀 여유는 없다. 그런데 “언제 도착하냐”만 알면 되는 경우가 의외로 많다.

아이코날 방정식(eikonal equation)은 매질의 국소 전파속도 $c(\mathbf{x})$ 가 주어졌을 때 파면의 도달 시간 $T(\mathbf{x})$ 가 만족하는 1계 비선형 편미분방정식이다.

$$\left|\nabla T(\mathbf{x})\right| = \frac{1}{c(\mathbf{x})} \equiv s(\mathbf{x}), \qquad T\big|_{\Gamma} = 0$$

$s$ 는 슬로우니스(slowness), $\Gamma$ 는 파원(점·곡선·곡면)이다. 이름은 그리스어 εἰκών(상, image)에서 왔고, 기하광학의 위상함수를 부르던 19세기 용어가 그대로 굳었다.

읽는 법은 간단하다. $T$ 의 등고선이 곧 파면이고, 등고선이 촘촘한 곳이 느린 곳이다. 진폭·간섭·회절 정보는 전부 버리고 “언제”만 남긴 것이라, 정보량을 극단적으로 줄인 대가로 헬름홀츠 방정식을 푸는 것보다 몇 자릿수 싸다. 파장이 매질 변화 척도보다 훨씬 짧을 때만 성립하는 근사지만, 지진 탐사·경로계획·레벨셋 방법의 재초기화처럼 그 조건이 맞아떨어지는 판이 넓다.

2. 파동방정식에서 나오는 길 — WKB 극한[편집]

출발점은 주파수영역 헬름홀츠 방정식 $\nabla^2 u + \omega^2 s^2 u = 0$ 이다. 여기에 고주파 전개(WKB 가설)

$$u(\mathbf{x}) = A(\mathbf{x})\,e^{i\omega T(\mathbf{x})}$$

를 넣고 $\omega$ 의 거듭제곱으로 정리하면 두 개의 방정식이 차례로 떨어진다.

$$O(\omega^{2}):\quad |\nabla T|^{2} = s^{2} \qquad\qquad O(\omega^{1}):\quad 2\,\nabla A\cdot\nabla T + A\,\nabla^{2}T = 0$$

앞쪽이 아이코날 방정식, 뒤쪽이 수송 방정식(transport equation)이다. 수송 방정식은 $\nabla\cdot(A^{2}\nabla T)=0$ 로 다시 쓸 수 있는데, 이는 광선 다발이 만드는 관을 따라 에너지 플럭스가 보존된다는 뜻이다. 즉 아이코날로 “언제”를, 수송으로 “얼마나 세게”를 얻는다.

버려진 항은 $\nabla^{2}A$ 다. 이 근사가 정당하려면 $|\nabla^{2}A| \ll \omega^{2}s^{2}A$, 실무적으로는 매질의 물성이 한 파장 안에서 거의 변하지 않아야 한다. 이 조건이 깨지는 곳 — 초점(caustic), 그림자 경계, 파장 규모의 산란체 — 에서는 진폭이 발산하거나 아이코날 해 자체가 물리와 어긋난다. 초점에서 $A\to\infty$ 가 되는 것은 기하광학의 고질병이고, 이를 손보려면 가우시안 빔이나 마슬로프 보정 같은 별도 장치가 필요하다.

3. 특성곡선 = 광선[편집]

$H(\mathbf{x},\mathbf{p}) = \tfrac12\left(|\mathbf{p}|^{2}-s^{2}\right)$, $\mathbf{p}=\nabla T$ 로 두면 아이코날은 $H=0$ 인 정상 해밀턴-야코비-벨만 방정식 계열이 되고, 특성곡선법의 정칙 절차가 그대로 돈다. 특성 방정식을 호길이 $\ell$ 로 다시 쓰면 익숙한 꼴이 나온다.

$$\frac{d}{d\ell}\!\left(s(\mathbf{x})\,\frac{d\mathbf{x}}{d\ell}\right) = \nabla s(\mathbf{x}), \qquad \frac{dT}{d\ell} = s(\mathbf{x})$$

이게 바로 광선 방정식이다. 아이코날 방정식의 특성곡선은 물리적인 광선(지진파라면 레이 경로)과 정확히 같고, 광선은 언제나 파면에 수직하게 진행한다. 레이 트레이싱이 광선을 하나씩 쏘는 라그랑주식 접근이라면, 아이코날을 격자에서 푸는 것은 같은 물리를 오일러식으로 푸는 것이다.

두 접근의 갈림길은 명확하다. 광선을 쏘면 임의 지점의 값을 얻으려고 조준(two-point ray tracing) 을 해야 하고, 광선이 성기게 지나가는 그림자 영역이 생긴다. 격자에서 풀면 도메인 전체에 빈틈없이 $T$ 가 채워지지만, 대신 다중 도달 시간을 표현할 수 없다 — 격자 한 점에 값은 하나뿐이기 때문이다.

4. 왜 점성해가 필요한가[편집]

문제는 광선이 서로 교차한다는 것이다. 속도장이 조금만 복잡해져도 서로 다른 경로를 거친 광선 두 개가 같은 점에 도달하고, 그러면 $T$ 는 다가 함수가 되려 한다. 격자 위의 단일값 함수로 강제하면 그 지점에서 $T$ 는 꺾이고($\nabla T$ 불연속), 고전적(모든 점에서 미분 가능한) 해는 존재하지 않는다.

더 나쁜 것은 약한 의미의 해가 무수히 많다는 사실이다. 1차원 $|T'|=1$, $T(0)=T(1)=0$ 만 봐도 톱니 모양으로 꺾인 립시츠 함수가 무한히 많고 전부 방정식을 거의 모든 점에서 만족한다. 어느 것이 “물리적인” 해인지 골라내는 기준이 필요하다.

그 기준이 크랜들과 라이온스(1983)의 점성해(viscosity solution)다. 원래 방정식에 인공 확산 $-\varepsilon\nabla^{2}T$ 를 넣어 매끈한 유일해를 얻은 뒤 $\varepsilon\to 0^{+}$ 극한을 취하는 것이 그 이름의 유래이고, 실제 정의는 시험함수와의 접촉 조건으로 확산 없이 서술한다. 부호 규약을 맞춰 잡으면 이 해는 최소 도달 시간, 즉 초동(first arrival)에 정확히 대응하고 유일하다. 아이코날 수치해법이 “안정적으로 수렴한다”고 말할 때의 수렴 대상은 언제나 이 점성해다.¹

여기서 반드시 기억할 대가가 하나 있다. 점성해는 초동만 준다. 저속도 층 아래를 돌아온 뒤늦은 도달파가 초동보다 훨씬 강한 진폭을 가지는 상황(지진 탐사에서 흔하다)에서도 아이코날 해는 약한 초동만 알려준다. 다중 도달 시간이 필요하면 파면 구성법(wavefront construction)이나 위상공간을 한 차원 늘려 레벨셋으로 다가 해를 추적하는 방법을 따로 써야 한다.

5. 상방향 이산화[편집]

점성해로 수렴하는 이산화의 핵심은 정보가 오는 쪽에서만 차분을 가져오는 것이다. 오셔·세시안(1988)의 상방향 도식, 루이·투랭(1992)이 점성해 수렴을 증명한 고두노프형 이산화가 표준이다. 2차원 균일격자 $h$ 에서

$$\max\!\left(D^{-x}T,\,-D^{+x}T,\,0\right)^{2} + \max\!\left(D^{-y}T,\,-D^{+y}T,\,0\right)^{2} = s_{ij}^{2}$$

이 식은 결국 격자점마다 작은 이차방정식 하나로 정리된다. 상류 이웃값 $T_1=\min(T_{i-1,j},T_{i+1,j})$, $T_2=\min(T_{i,j-1},T_{i,j+1})$ 를 잡고

$$(T-T_1)^{2} + (T-T_2)^{2} = h^{2}s^{2} \;\Longrightarrow\; T = \frac{T_1+T_2+\sqrt{2h^{2}s^{2}-(T_1-T_2)^{2}}}{2}$$

큰 근을 취한다. 판별식이 음수면(=두 이웃의 차이가 $hs$ 보다 커서 한 방향에서만 정보가 오는 경우) 1차원 갱신 $T=\min(T_1,T_2)+hs$ 로 후퇴한다. 이 큰 근 · 인과성 검사 두 줄이 아이코날 솔버의 심장이고, 나머지는 전부 “이 갱신을 어떤 순서로 도느냐”의 문제다.

6. 고속행진법[편집]

세시안(1996)의 고속행진법(Fast Marching Method, FMM)은 순서 문제를 다익스트라 알고리즘처럼 푼다. 상방향 갱신식은 자기보다 작은 값에만 의존하므로, 값이 작은 것부터 확정해 나가면 확정된 값은 두 번 다시 바뀌지 않는다.

• 격자점을 확정(Accepted) · 후보(Trial) · 미방문(Far) 세 상태로 나눈다.
• 후보 집합을 최소 힙에 넣고, 가장 작은 후보를 꺼내 확정한다.
• 확정된 점의 이웃을 위 이차방정식으로 갱신하고 힙을 정리한다.

한 점이 정확히 한 번 확정되므로 전체는 한 번의 통과로 끝나고, 비용은 힙 연산이 지배해 $O(N\log N)$ 이다. 격자 크기와 무관하게 반복이 없다는 점이 강력하다. 최적제어 쪽에서 치치클리스(1995)가 사실상 같은 구조를 독립적으로 얻었기 때문에, 문헌에 따라 두 사람을 함께 인용한다.

FMM에도 함정은 있다.

• 비정렬 격자에서 인과성이 깨진다. 삼각형이 둔각이면 갱신이 “더 큰 값”에 의존하게 되어 단조성이 무너진다. 킴멜·세시안(1998)의 펼침(unfolding)/가상 간선 분할이 표준 대응책이다.
• 비등방 문제로 확장되지 않는다. $H$ 가 방향에 의존하면 최소값 우선 순서가 더 이상 정답이 아니다. 세시안·블라디미르스키의 정렬 상방향 방법(ordered upwind method)은 비등방 비율에 비례해 스텐실을 넓혀 이 문제를 푼다. 스텐실 폭이 비등방성에 비례해 커지므로 공짜는 아니다.
• 병렬화가 어렵다. 힙은 본질적으로 순차적이다. GPU에서는 힙을 버리고 활성 리스트를 쓰는 고속 반복법(FIM) 계열로 갈아탄다.

7. 고속 스위핑법[편집]

고속 스위핑법(Fast Sweeping Method, FSM)은 정반대의 발상이다. 정렬을 포기하고 가우스-자이델식으로 격자를 훑되, 훑는 방향을 번갈아 바꾼다. 2차원이면 $(i{\uparrow}j{\uparrow})$, $(i{\downarrow}j{\uparrow})$, $(i{\uparrow}j{\downarrow})$, $(i{\downarrow}j{\downarrow})$ 네 방향, $d$ 차원이면 $2^{d}$ 방향이다. 각 스윕이 특성곡선 방향의 한 사분면(팔분면)을 담당하고, 갱신은 항상 $T \leftarrow \min(T,\,\text{새 값})$ 으로 단조 감소시킨다.

자오(2005)가 증명한 핵심은 필요한 스윕 횟수가 격자 간격 $h$ 에 무관한 상수라는 것이다. 따라서 총 비용은 $O(N)$ 으로, 정렬도 힙도 없다. 단순 루프라 캐시 친화적이고 상수도 작아, 속도가 매끈한 문제에서는 FMM보다 실측이 몇 배 빠른 경우가 흔하다.

물론 그 “상수”는 속도장이 얼마나 광선을 휘게 하느냐에 붙어 있다. 특성곡선이 도메인 안에서 방향을 여러 번 바꾸는 복잡한 매질에서는 스윕 횟수가 늘어난다. 요약하면 이렇다.

항목	고속행진법	고속 스위핑법
복잡도	$O(N\log N)$, 매질 무관	$O(N)$, 상수가 매질 의존
자료구조	최소 힙 필요	배열 루프뿐
최악의 경우	예측 가능	굴곡 심한 속도장에서 스윕 증가
병렬화	어려움 (힙이 순차적)	상대적으로 쉬움 (영역 분할·대각선 스윕)

실무에서는 둘 다 구현해 두고 문제에 따라 고르는 것이 정답에 가깝다.

8. 정확도의 함정 — 파원 특이성[편집]

여기서 초심자가 반드시 한 번 데인다. 1차 상방향 도식은 $O(h)$ 인데, 2차 도식으로 바꿔도 전체 오차가 1차로 주저앉는 일이 생긴다. 원인은 파원이다.

점 파원 근처에서 해는 $T \approx s_0|\mathbf{x}-\mathbf{x}_0|$ 이라 값 자체는 멀쩡하지만 2계 도함수가 발산한다. 유한차분의 절단오차는 고계 도함수에 비례하므로 파원 주변 몇 셀에서 국소 오차가 폭발하고, 그 오차가 특성곡선을 타고 도메인 전체로 실려 나간다. 대응책은 두 갈래다.

• 파원 주변 작은 공을 해석해로 직접 채워 넣고 그 바깥에서만 차분을 시작한다. 간단하지만 속도가 변하는 매질에서는 채우는 반경 선정이 애매하다.
• 인수분해 아이코날(factored eikonal): $T = T_{0}\,\tau$ 로 두고 $T_{0}$ 를 균질 매질의 알려진 해로 잡으면, 특이성이 $T_0$ 에 흡수되어 남는 $\tau$ 는 파원에서 매끄럽다. $\tau$ 에 대한 방정식을 FMM/FSM으로 풀면 고차 도식이 제 차수를 낸다. 오늘날 지진 탐사 쪽 실무 표준에 가깝다.

또 하나. 격자 수렴 검사를 할 때 오차 노름을 무엇으로 잴지 주의해야 한다. 점성해는 꺾인 점(shock)을 가지므로 $L^\infty$ 오차는 그 근방에서 잘 떨어지지 않는다. $L^1$ 로 재면 깔끔하게 수렴 차수가 보인다. 수치 소산과 분산에서 다루는 이야기가 여기서도 그대로 반복된다.

9. 응용[편집]

초동주시 계산과 지진 토모그래피. 지하 속도 모델에서 각 발파점의 초동 도달 시간을 계산하는 것이 아이코날의 최대 고객이다. 비데일(1988)의 확장 사각형 유한차분법이 시작이었으나 속도 대비가 큰 곳에서 인과성이 깨져 불안정했고, 이를 FMM이 정리했다. 관측 주시와 계산 주시의 차이를 속도 모델의 갱신으로 되돌리는 것이 지진파 토모그래피이며, 감도 커널이 광선을 따라 놓이므로 전형적인 부적절(ill-posed) 역문제가 된다. 초동만 쓴다는 제약은 위에서 말한 그대로다.

경로계획. 속도 $c(\mathbf{x})$ 를 “그 지점을 지나기 쉬운 정도”로 해석하면 $T$ 는 비용-투-고(cost-to-go)가 되고, $-\nabla T$ 를 따라 내려가면 최적 경로가 나온다. 격자 위 A*가 8방향 연결 때문에 진행 방향을 45° 배수로 양자화하는 격자화 오차를 겪는 반면, 아이코날 해는 연속 방향을 제대로 낸다. 내비게이션 메시의 흐름장, 군중 시뮬레이션의 연속체 모형, 최적 제어의 값함수가 모두 같은 방정식의 다른 이름이다.

거리함수와 재초기화. $c\equiv 1$ 이면 $|\nabla T|=1$ 이라 $T$ 는 곧 파원까지의 거리다. 레벨셋 방법에서 이류로 망가진 $\phi$ 를 부호거리함수로 되돌리는 재초기화를, 가상시간 PDE 대신 0-등고선에서 FMM을 한 번 행진시켜 해결하는 것이 널리 쓰인다. 계면 위치를 흔들지 않아 질량 누출이 덜하다는 것이 장점. 다만 순수 유클리드 거리만 필요하다면 거리 변환의 정확 알고리즘이 더 싸고 정확하다 — 아이코날은 속도가 변할 때 값을 하는 도구다. 수학적 형태학의 측지 거리 연산도 같은 자리에 놓인다.

그 밖에. 삼각망 위 측지선 계산, 음영으로부터 형상 복원(shape-from-shading), 반도체 공정의 식각·증착 표면 진화, 영상 분할의 측지 활성 윤곽, 마칭 큐브로 뽑을 등가면의 사전 계산까지 — “가장 빠른 도달”이라는 말이 나오면 대개 이 방정식이 뒤에 있다.²

10. 여담[편집]

• 이름이 어려워 보이지만 하는 일은 “번짐”이다. 실제로 FMM을 처음 구현한 사람들의 소감은 대체로 “다익스트라를 격자에서 돌린 것 아니냐”인데, 절반은 맞다. 나머지 절반, 즉 이차방정식 갱신이 격자화 오차를 없애는 부분이 이 방법의 알맹이다.
• 속도장에 0이나 음수가 들어가면 즉시 터진다. 슬로우니스가 무한대가 되기 때문. 장애물은 속도 0이 아니라 매우 작은 양수로 넣거나 아예 도메인에서 빼는 것이 실무 관행이다.³
• 아이코날은 충격파 문제와 사촌이다. 보존형 방정식의 충격이 도달 시간의 꺾임에 대응하고, 전변분 감소·WENO 도식 계열의 고차 상방향 기법이 그대로 이식된다. 아이코날용 고차 도식 논문의 참고문헌이 CFD 논문과 겹치는 것은 우연이 아니다.⁴

11. 관련 문서[편집]

• 레벨셋 방법 · 거리 변환 · 수학적 형태학
• 특성곡선법 · 해밀턴-야코비-벨만 방정식 · 최적 제어
• 헬름홀츠 방정식 · 레이 트레이싱 · 물리광학
• WENO 도식 · ENO 도식 · 수치 소산과 분산 · 유한차분법
• 역문제 · 동적 계획법 · 내비게이션 메시 · 군중 시뮬레이션
• 변환광학 · 등각사상

12. Footnotes[편집]

1. 이름 때문에 “점성 유체와 관계가 있나” 하는 오해가 종종 생기는데 전혀 없다. 증명 기법에서 인공 점성을 썼다는 역사적 흔적일 뿐이고, 정의 자체에는 $\varepsilon$ 이 나오지 않는다. 참고로 크랜들과 라이온스의 1983년 논문은 그 뒤 20년간 이 바닥 인용 랭킹을 독식했다. ↩

2. 세시안 본인이 이 방법으로 반도체 공정 시뮬레이션 회사를 차렸다. 수치해석 논문이 그대로 제품이 되는 흔치 않은 사례이고, 그래서 FMM 관련 문헌에는 유독 특허 각주가 많다. ↩

3. 진짜 0을 넣으면 힙에 무한대가 들어가고, 부동소수점 비교가 어긋나면서 힙 불변식이 깨져 프로그램이 조용히 틀린 답을 낸다. 발산하지 않고 그럴듯하게 틀린 답을 내는 것이 이 버그의 가장 악질적인 부분이다. ↩

4. 반대 방향의 이식도 있다. 정상 상태 아이코날을 잘 푸는 스위핑 아이디어가 정상 상태 이류-확산 문제의 전처리기로 되돌아간 사례가 여럿이다. 결국 다들 상류가 어디냐를 놓고 싸우는 중이다. ↩