A Tutorial on the Analytical Solution of Lamb’s Problem and Its Numerical Implementation

Sangin Cho; Sukjoon Pyun

doi:10.7582/GGE.2026.29.2.131

Preview

Geophysics and Geophysical Exploration. 31 May 2026. 131-142
https://doi.org/10.7582/GGE.2026.29.2.131

A Tutorial on the Analytical Solution of Lamb’s Problem and Its Numerical Implementation

Lamb 문제 해석해와 수치 구현에 대한 해설

Sangin Cho¹

Sukjoon Pyun¹^*

조 상인¹

편 석준¹^*

¹Department of Energy Resources Engineering, Inha University

¹인하대학교 에너지자원공학과

^{*Corresponding Author}

License (open-access, http://creativecommons.org/licenses/by-nc/4.0/):

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

ABSTRACT

Lamb’s problem is a core fundamental topic for understanding wave phenomena in seismology and geophysical exploration. This paper details the derivation of the analytical solution and its numerical implementation for a 2D semi-infinite medium subjected to a vertical point load on the surface. We derive the solution in the frequency-wavenumber domain using Helmholtz decomposition and Fourier transforms and present a step-by-step algorithm to reconstruct the time-space wavefield via the discrete Fast Fourier Transform (FFT). Comparing the implemented analytical wavefields with finite element method (FEM) modeling results confirmed that the travel times and waveforms of major seismic events match with high accuracy. This paper will serve as a useful guide for beginners to practically understand the basic concepts of elastic analytical solutions and their numerical implementation.

Keywords

Lamb’s Problem

Analytical Solution

Elastic Wave Equation

Numerical Implementation

Finite Element Method

Lamb 문제는 지진학 및 물리탐사 분야에서 파동 현상을 이해하기 위한 핵심적인 기초 주제이다. 이 해설 논문은 2차원 반무한 매질 표면에 수직 점하중이 작용하는 경우에 대한 해석해 유도와 수치적 구현 과정을 상세히 다룬다. 헬름홀츠 분해와 푸리에 변환을 이용해 주파수-파수 영역의 해를 유도하고, 이산 고속 푸리에 변환(FFT)을 통해 시간-공간 영역 파동장으로 복원하는 알고리즘을 단계별로 제시한다. 구현된 해석해 파동장을 유한요소법(FEM) 모델링 결과와 비교한 결과, 주요 파동의 도달 시간과 파형이 매우 정확히 일치함을 확인했다. 이 논문은 입문자들이 탄성파 해석해의 기본 개념과 수치적 구현 과정을 실무적으로 이해하는 데 유용한 지침이 될 것이다.

키워드

Lamb 문제

해석해

탄성파동방정식

수치구현

유한요소법

MAIN

서 론
이 론
수치 예제
결 론

서 론

탄성 매질에서의 파동 전파를 해석적으로 표현하려는 연구는 이론 지진학과 물리탐사 분야에서 오랫동안 중요한 주제로 다루어져 왔다. 물리적 현상을 수학적으로 명확히 표현하는 해석해는 파동 전파 현상을 이해하는 데 중요한 역할을 하며, 수치적 기법으로 구한 해와 달리 문제의 매개변수가 결과에 미치는 영향을 명확히 파악할 수 있다. 해석해는 탄성파 자료의 모델링 및 해석(Aki and Richards, 2002), 표면파와 실체파의 전파 특성 분석(Graff, 1975), 송신원과 수신기 배치에 따른 파동장 거동 이해(Červený, 2001; Hadley and Helmberger, 1977), 수치 모델링 및 역산 기법의 검증(Kausel, 2013)뿐만 아니라 지반 응답 해석(site response analysis) 및 부지 증폭 효과 평가(Bard and Bouchon, 1985) 등 다양한 분야에서 중요한 역할을 한다. 이와 같은 이유로 기초적인 탄성파동 해석해에 대한 이해는 탄성파 관련 연구를 수행하는 데 중요한 기반이 된다. 본 해설 논문에서는 Lamb 문제(Lamb’s problem)를 중심으로 탄성파 해석해의 유도 과정과 수치 구현 방법을 설명한다.

Lamb 문제는 등방성(isotropic)이고 균질(homogeneous)한 반무한 탄성 매질에서 동적 점 하중(dynamic point load)에 의해 발생하는 파동장을 해석적으로 구하는 문제이다. 이 문제의 기원은 19세기 후반으로 거슬러 올라간다. Rayleigh (1885)는 탄성 반무한 공간에서 파동 전파 문제의 새로운 해를 제시하며, 고체 표면을 따라 전파하는 표면파의 존재를 규명하였다. 이를 통해 탄성체의 무응력(traction-free) 경계 조건에서 표면파의 운동 특성이 기존에 알려진 실체파(body wave)와 다름을 입증하였다. 이후 Lamb (1904)은 지진 현상이 지구 표면에서 나타나는 방식을 설명하기 위해 외부 힘이 작용할 때 자유면의 운동, 즉 경계에서의 파동방정식의 해를 구하고자 하였다.

Lamb의 연구는 반무한 탄성 매질에서 선 하중(line load)에 의해 발생하는 2차원 파동장을 해석적으로 구하는 문제로 시작된다. 그는 헬름홀츠 분해(Helmholtz decomposition)와 푸리에 적분 변환(Fourier integral transform)을 이용하여 2차원 단일 주파수 해(monochromatic solutions)를 유도하였다. 이 해는 적분 경로 상의 특이점을 포함하는 이상 적분(improper integral) 형태로 표현된다. 이를 3차원으로 확장하여 조화(time-harmonic) 파동장 해를 원통 좌표계(cylindrical coordinates)에서 유도하였다. 또한 시간영역 해(transient solutions)를 구하기 위해 푸리에 합성(Fourier synthesis)을 적용하였으며, 이 과정에서 각주파수에 대한 추가적인 푸리에 적분이 수행되어 최종 해는 이중 적분(double integrals) 형태로 표현되었다(Emami and Eskandari-Ghadi, 2019). Lamb이 제시한 문제는 외력이 수직 방향으로 작용할 때 매질 경계에서의 운동을 설명하지만, 수평 또는 기울어진 힘에 의한 파동 전파나 매질 내부의 특정 지점에서의 해는 다루지 않았다. 이에 따라 보다 일반화된 해석해를 얻기 위한 연구가 활발히 진행되었다(Miller and Pursey, 1954; Pekeris, 1955a; De Hoop, 1960). 또한, 시간영역 파동장을 얻기 위해 주파수 영역 해를 푸리에 합성하는 과정은 수학적으로 복잡하고 계산 효율성이 낮다는 한계가 있어, 이를 개선하기 위한 다양한 복소 해석적 기법이 제안되었다(Cagniard, 1962; De Hoop, 1960).

Lamb 문제는 송신원과 수신기의 위치에 따라 유형 1, 유형 2, 유형 3으로 구분된다(Hu and Jiang, 2024). 유형 1은 송신원과 수신기가 모두 지표에 존재하는 경우, 유형 2는 송신원이 매질 내부에 있고 수신기가 지표에 존재하는 경우, 유형 3은 송신원과 수신기 모두 매질 내부에 존재하는 경우를 나타낸다. 유형 1에 대해 Pekeris (1955a)와 Chao (1960)는 특정 포아송 비( $ν = 0.25$ )를 갖는 매질에서 각각 수직 및 수평 점 송신원에 대한 닫힌 해(closed-form solution)를 유도하였다. 이후 Mooney (1974)는 이를 임의의 포아송 비를 가지는 매질로 일반화하여 수직 방향 외력에 대한 해를 적분 형태로 제시하였으며, Richards (1979)와 Kausel (2013)은 유형 1에 대해 완전한 닫힌 해를 제시하였다.

유형 2와 유형 3의 경우에는 보다 복잡한 경계 조건과 내부 송신원 효과를 다루기 위해 다양한 해석 기법이 개발되었다. Lapwood (1949), Pinney (1954), Pekeris (1955b) 등의 연구가 이어졌으며, 특히 Cagniard (1962)가 제안하고 De Hoop (1960)이 발전시킨 Cagniard–de Hoop 방법은 Lamb 문제를 해석하는 대표적인 접근법으로 자리 잡았다. 이 방법은 라플라스 변환과 복소 해석학적 기법을 이용하여 파동 전파 문제의 적분 해를 효율적으로 계산할 수 있도록 하며, 이후 다양한 Lamb 문제의 해석해 유도에 널리 활용되었다. Johnson (1974)은 Cagniard–de Hoop 방법을 이용하여 반무한 공간 내부 및 경계에서 수직·수평 점 하중이 작용하는 경우에 대한 적분 형태의 해를 제시하였다. 이후 Feng and Zhang (2018, 2020a)은 유형 2와 유형 3에 대한 닫힌 형태의 해를 유도하고, 이동하는 송신원에 대한 정확한 닫힌 해를 제시하였다(Feng and Zhang, 2020b). 그러나 이들의 해는 포아송 비가 특정한 범위 ( $0 < ν < 0.2631$ )에 존재할 때만 유효하다는 한계를 가진다. 이후 Hu and Jiang (2024)는 이를 확장하여 임의의 포아송 비에 대해 적용 가능한 닫힌 해를 제시하였다.

이 해설 논문에서는 앞서 언급한 여러 풀이 방법 중 Lamb (1904)이 제시한 고전적인 적분 표현 기반 접근법을 중심으로 해석해를 구하는 과정을 살펴보고 이를 코드로 구현하는 방법을 설명한다. 특히 기존 연구들이 해석해의 수학적 유도에 초점을 맞춘 것과 달리, 본 논문에서는 실제 코드 구현과 시간–공간 영역 파동장 계산 과정에 중점을 둔다. 최근 컴퓨터 프로그래밍 환경의 발전으로 푸리에 변환 기반 계산은 비교적 효율적으로 수행될 수 있으며, 이를 통해 해석해 기반 파동장 계산의 실제 구현 가능성이 크게 향상되었다. 이 연구에서는 2차원 공간에서 유형 1에 해당하는 Lamb 문제에 대한 적분 형태의 해를 유도하고, 이를 기반으로 시간–공간 영역의 파동장을 계산하였다. 이때 자유표면에는 수직 방향의 점하중(vertical point load)이 작용한다고 가정하였으며, 지표에 도달하는 파동장은 공간 푸리에 변환을 통해 평면파 성분의 중첩(superposition of plane waves) 형태로 표현된다. 구현된 해석해는 수치 모델링 결과와 비교하여 타당성을 검증하고자 한다.

이 론

우선 지배 방정식인 힘 평형 방정식(Force balance equation)에서 체적력(body force)이 없다고 가정하고, 후크의 법칙(Hooke’s law)을 적용하여 탄성 파동 방정식을 도출한다. 2차원 탄성파 문제를 고려할 때 파동장의 두 성분을 각각 $u$ 와 $w$ 라 하면, 지배 방정식은 다음과 같다.

(1)

ρ \frac{\partial^{2} u}{\partial t^{2}} = (λ + μ) \frac{\partial ∆}{\partial x} + μ \nabla^{2} u

(2)

ρ \frac{\partial^{2} w}{\partial t^{2}} = (λ + μ) \frac{\partial ∆}{\partial z} + μ \nabla^{2} w

여기서 $λ$ 와 $μ$ 는 Lamé 상수이며, $ρ$ 는 밀도를 의미한다. 헬름홀츠 분해(Helmholtz decomposition)는 파동장을 스칼라 포텐셜과 벡터 포텐셜로 분리하여, 각각의 성분이 독립적으로 해석될 수 있도록 한다. 헬름홀츠 포텐셜을 이용하면 파동장 벡터 $u$ 는 다음과 같다.

(3)

u = \nabla ϕ + \nabla \times ψ

여기서 $ϕ$ 는 압축파(compressional wave) 성분에 대응하는 스칼라 포텐셜이며, $ψ$ 는 전단파(shear wave) 성분을 나타내는 벡터 포텐셜 함수이다. 헬름홀츠 포텐셜로부터 2차원 파동장의 두 성분 $u$ (수평 성분)와 $w$ (수직 성분)는 다음과 같이 표현된다.

(4)

u = \frac{\partial ϕ}{\partial x} + \frac{\partial ψ}{\partial z}

(5)

w = \frac{\partial ϕ}{\partial z} - \frac{\partial ψ}{\partial x}

식 (4)와 (5)에서 정의한 파동장과 포텐셜 함수 사이의 관계를 이용하면 파동방정식과 경계조건을 포텐셜 함수로 표현할 수 있다.

(6)

\nabla^{2} ϕ = \frac{1}{α^{2}} \frac{\partial^{2} ϕ}{\partial t^{2}}

(7)

\nabla^{2} ψ = \frac{1}{β^{2}} \frac{\partial^{2} ψ}{\partial t^{2}}

(8)

P_{z x} = μ [\frac{\partial w}{\partial x} + \frac{\partial u}{\partial z}] = μ [2 \frac{\partial^{2} ϕ}{\partial x \partial z} - \frac{\partial^{2} ψ}{\partial x^{2}} + \frac{\partial^{2} ψ}{\partial z^{2}}]

(9)

P_{z z} = λ [\frac{\partial u}{\partial x} + \frac{\partial w}{\partial z}] + 2 μ \frac{\partial w}{\partial z} = λ [\frac{\partial^{2} ϕ}{\partial x^{2}} + \frac{\partial^{2} ϕ}{\partial z^{2}}] + 2 μ [\frac{\partial^{2} ϕ}{\partial z^{2}} - \frac{\partial^{2} ψ}{\partial x \partial z}]

여기서 $α$ 와 $β$ 는 각각 P파와 S파 속도를 의미하며, $P_{z x}$ 와 $P_{z z}$ 는 각각 전단 응력과 수직 응력 성분을 나타낸다. 식 (6)과 (7)을 풀기 위해, 포텐셜 함수 $ϕ$ , $ψ$ 를 시간조화 평면파(time-harmonic plane wave) 함수로 가정할 수 있다.

(10)

ϕ = A e^{- i (k_{x} x + k_{z P} z)} e^{i ω t}

(11)

ψ = B e^{- i (k_{x} x + k_{z S} z)} e^{i ω t}

여기서 $ω$ 는 각주파수이며, $k_{x}$ 는 $x$ 방향에 대한 수평 파수, $k_{z P}$ 와 $k_{z S}$ 는 각각 P파와 S파 포텐셜에 대한 $z$ 방향 파수이다. 반무한 탄성매질은 $x$ 방향으로 균질하므로, 파동장은 $x$ 방향에 대해 $e^{- i k_{x} x}$ 형태의 푸리에 성분으로 표현할 수 있다. 따라서 P파와 S파 모두 동일한 수평 방향 파수 $k_{x}$ 를 공유한다. 반면, 두 파는 서로 다른 전파 속도를 가지므로 전체 파수 크기가 다르며, 이에 따라 수직 방향 파수 성분은 서로 다르게 결정된다.

(12)

k_{P} = \frac{ω}{α}

(13)

k_{S} = \frac{ω}{β}

P파와 S파의 파수 $k_{P}$ , $k_{S}$ 를 식 (12), (13)과 같이 나타내면, 수직 방향 파수 성분 $k_{z P}$ 와 $k_{z S}$ 는 각각 다음과 같은 관계를 만족한다.

(14)

k_{z P} = \sqrt{k_{P}^{2} - k_{x}^{2}} = k_{x} \sqrt{\frac{c^{2}}{α^{2}} - 1}

(15)

k_{z S} = \sqrt{k_{S}^{2} - k_{x}^{2}} = k_{x} \sqrt{\frac{c^{2}}{β^{2}} - 1}

여기서 $c = \frac{ω}{k_{x}}$ 는 $x$ 방향 위상속도이다. 식 (14), (15)의 수직 방향 파수 성분을 식 (10), (11)에 대입하여 다음과 같이 정리한다.

(16)

ϕ = A e^{- i k_{x} (x + \sqrt{\frac{c^{2}}{α^{2}} - 1} z)} e^{i ω t} = A e^{- i k_{x} x - p z} e^{i ω t}

(17)

ψ = B e^{- i k_{x} (x + \sqrt{\frac{c^{2}}{β^{2}} - 1} z)} e^{i ω t} = B e^{- i k_{x} x - q z} e^{i ω t}

여기서 $p (= k_{x} \sqrt{\frac{c^{2}}{α^{2}} - 1})$ 와 $q (= k_{x} \sqrt{\frac{c^{2}}{β^{2}} - 1})$ 는 각각 P파와 S파에 대한 수직 방향 감쇠 계수를 의미한다.

지표면에서 수직 방향으로 점하중이 작용하는 경우, 이는 $z = 0$ 에서 수직 응력 형태의 점하중 경계 조건으로 표현할 수 있다(식 (18)).

(18)

P_{z z} = f (t) δ (x)

여기서 $f (t)$ 는 시간에 따른 송신원 함수(source function)이며, $δ (\cdot)$ 는 디락–델타(dirac-delta) 함수이다. 수평 방향으로 작용하는 외력은 없다고 가정하였으므로 $z = 0$ 에서 다음 조건 또한 만족해야 한다.

(19)

P_{z x} = 0

식 (16), (17)을 식 (18), (19)의 경계조건에 대입하면 계수 $A$ , $B$ 에 대한 연립방정식을 얻을 수 있다.

(20)

(2 {k_{x}}^{2} - {k_{x}}^{2} \frac{c^{2}}{β^{2}}) A - 2 i k_{x} q B = \frac{1}{μ}

(21)

2 i k_{x} p A + (2 {k_{x}}^{2} - {k_{x}}^{2} \frac{c^{2}}{β^{2}}) B = 0

식 (20), (21)의 방정식을 연립하여 풀면, 미지수 $A$ , $B$ 는 다음과 같이 결정된다.

(22)

A = \frac{1}{μ} \frac{2 {k_{x}}^{2} - {k_{x}}^{2} \frac{c^{2}}{β^{2}}}{{(2 {k_{x}}^{2} - {k_{x}}^{2} \frac{c^{2}}{β^{2}})}^{2} - 4 {k_{x}}^{2} p q}

(23)

B = - \frac{1}{μ} \frac{2 i k_{x} p}{{(2 {k_{x}}^{2} - {k_{x}}^{2} \frac{c^{2}}{β^{2}})}^{2} - 4 {k_{x}}^{2} p q}

이를 적용하여 포텐셜 함수 $ϕ (x, t)$ 와 $ψ (x, t)$ 를 다음과 같이 적분식으로 나타낼 수 있다.

(24)

ϕ (x, t) = \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} \frac{1}{μ} \frac{2 {k_{x}}^{2} - {k_{x}}^{2} \frac{c^{2}}{β^{2}}}{{(2 {k_{x}}^{2} - {k_{x}}^{2} \frac{c^{2}}{β^{2}})}^{2} - 4 {k_{x}}^{2} p q} e^{- i k_{x} x} e^{i ω t} d k_{x} d ω

(25)

ψ (x, t) = \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} - \frac{1}{μ} \frac{2 i k_{x} p}{{(2 {k_{x}}^{2} - {k_{x}}^{2} \frac{c^{2}}{β^{2}})}^{2} - 4 {k_{x}}^{2} p q} e^{- i k_{x} x} e^{i ω t} d k_{x} d ω

식 (24), (25)의 포텐셜 함수를 식 (4), (5)에 대입하여 파동장의 두 성분을 다음과 같이 표현할 수 있다.

(26)

u = \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} - \frac{1}{μ} \frac{k_{x} (2 {k_{x}}^{2} - {k_{x}}^{2} \frac{c^{2}}{β^{2}} - 2 p q)}{{(2 {k_{x}}^{2} - {k_{x}}^{2} \frac{c^{2}}{β^{2}})}^{2} - 4 {k_{x}}^{2} p q} e^{- i k_{x} x} e^{i ω t} d k_{x} d ω

(27)

w = \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} \frac{1}{μ} \frac{p {k_{x}}^{2} \frac{c^{2}}{β^{2}}}{{(2 {k_{x}}^{2} - {k_{x}}^{2} \frac{c^{2}}{β^{2}})}^{2} - 4 {k_{x}}^{2} p q} e^{- i k_{x} x} e^{i ω t} d k_{x} d ω

식 (26)과 (27)에서 나타난 바와 같이 파동방정식의 해는 주파수와 파수 영역에서 정의된 함수의 푸리에 변환 형태임을 알 수 있다. 이 커널 함수를 주파수–파수 영역에서 계산하고 푸리에 역변환함으로써 시간–공간 영역의 파동장을 복원할 수 있다.

한편, 점하중이 작용하는 $x$ 축 상의 위치가 원점이 아닌 경우를 고려하기 위해 경계조건을 다음 식과 같은 형태로 확장하였다.

(28)

P_{z z} = f (t) δ (x - x_{0})

이는 지표면( $z = 0$ )에서 특정한 위치( $x = x_{0}$ )에 점 하중이 작용하는 경우를 의미한다. 따라서 기존 경계조건으로부터 얻어진 식 (20)은 다음과 같이 수정되어야 한다.

(29)

(2 {k_{x}}^{2} - {k_{x}}^{2} \frac{c^{2}}{β^{2}}) A - 2 i k_{x} q B = \frac{e^{i k_{x} x_{0}}}{μ}

식 (29)를 이용하면 변위에 대한 최종 해는 다음과 같이 표현된다.

(30)

u = \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} - \frac{i e^{i k_{x} x_{0}}}{μ} \frac{k (2 {k_{x}}^{2} - {k_{x}}^{2} \frac{c^{2}}{β^{2}} - 2 p q)}{{(2 {k_{x}}^{2} - {k_{x}}^{2} \frac{c^{2}}{β^{2}})}^{2} - 4 {k_{x}}^{2} p q} e^{- i k_{x} x} e^{i ω t} d k_{x} d ω

(31)

w = \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} \frac{e^{i k_{x} x_{0}}}{μ} \frac{p {k_{x}}^{2} \frac{c^{2}}{β^{2}}}{{(2 {k_{x}}^{2} - {k_{x}}^{2} \frac{c^{2}}{β^{2}})}^{2} - 4 {k_{x}}^{2} p q} e^{- i k_{x} x} e^{i ω t} d k_{x} d ω

여기서 점하중의 위치 변화는 푸리에 영역(주파수–파수 영역)에서 위상 항 $e^{i k_{x} x_{0}}$ 를 곱하는 형태로 구현된다.

수치 예제

매질의 경계(지표면)에서 수직 방향의 점하중이 작용하는 경우 발생하는 파동장의 수평 및 수직 변위 성분은 식 (30)과 (31)을 수치적으로 계산하여 얻을 수 있다. 특히, 포트란(Fortran), C, 파이썬(Python) 등의 프로그래밍 언어에서는 푸리에 변환 연산이 함수 또는 라이브러리 형태로 제공되므로 해석해를 비교적 용이하게 계산할 수 있다.

먼저 모델링에 필요한 파라미터를 설정한다. 매질의 특성을 정의하기 위해 필요한 파라미터는 매질의 P파 속도( $α$ ), S파 속도( $β$ ), 밀도( $ρ$ )이다. 또한, 공간 격자 간격( $∆ x$ ), 공간 격자 개수( $n_{x}$ ), 최대 주파수( $f_{m a x}$ ), 관측 시간( $t_{m a x}$ ), 시간 영역 출력 파형의 시간 간격( $∆ t_{r}$ ), 점하중의 수평 위치( $x_{0}$ ) 등을 정의해야 한다. 제로 패딩(zero padding)은 주파수 영역 배열의 끝에 0을 추가하여 주파수 축 샘플 개수를 증가시키는 방법으로, 시간 영역 신호의 샘플 간격을 감소시켜 보간 효과를 제공한다. 추가로 푸리에 변환 계산에서 발생하는 두루마리 왜곡(wrap-around artifact)을 완화하기 위한 감쇠 계수(damping coefficient)를 설정해야 한다. 두루마리 왜곡은 푸리에 변환이 입력 신호를 주기 함수(periodic function)로 가정하기 때문에 발생하는 비물리적 현상으로, 계산 영역의 한쪽 끝을 벗어난 파동이 반대쪽 끝으로 다시 나타나는 현상을 의미한다. 본 논문에서는 이를 완화하기 위해 복소 주파수(complex frequency)를 적용하였으며, 복소 각주파수를 $Ω = ω - i γ$ 로 정의하였다.

해석해 계산에 필요한 추가 파라미터들은 입력 파라미터를 이용하여 계산되며, 그 정의는 Table 1에 정리하였다. 프로그래밍 언어로 계산을 수행할 때는 연속적인 값을 사용하는 것이 불가능하기 때문에 이산화된 파라미터를 사용한다. 특히 영역 변환 시에도 신호가 훼손되지 않도록 샘플링 이론에 따라 시간 간격( $∆ t$ )과 공간 간격( $∆ x$ )을 결정해야 한다.

Table 1

Definitions and corresponding equations for the principal parameters used in the analytical solution calculation.

Parameter	Formula
Frequency interval ( $∆ f$ ) [Hz]	$\frac{1}{t_{m a x}}$
Number of frequency samples ( $n_{f}$ )	$i n t (\frac{f_{m a x}}{∆_{f}})$
Time interval ( $∆_{t}$ ) [sec]	$\frac{1}{2 * f_{m a x}}$
Number of time samples ( $n_{t}$ )	$2 * (n_{f} - 1)$
Horizontal model length ( $x_{m a x}$ ) [m]	$(n_{x} - 1) * ∆ x$
Wavenumber interval ( $∆_{k}$ ) [m^-1]	$\frac{1}{x_{m a x}}$
Maximum wavenumber ( $k_{m a x}$ ) [m^-1]	$\frac{1}{2 * ∆ x}$
Number of wavenumber samples ( $n_{k}$ )	$i n t (\frac{k_{m a x}}{∆ k})$
Number of interpolated time samples ( $n_{r t}$ )	$i n t (\frac{t_{m a x}}{∆ t_{r}})$

샘플 개수는 반드시 0 또는 양의 정수 값이므로 최대 크기와 간격으로 나눈 후 정수화 하는 과정이 필수적이다. 이때 프로그래밍 언어에 따라 정수화 방식이 다르므로 파라미터 설정에 유의해야 한다.

본격적인 파동장 계산에 앞서 송신파형을 결정한다. 여기서는 탄성파 모델링에서 주로 사용되는 리커 파형(Ricker wavelet)을 사용하였다. 최대 주파수의 1/3을 피크 주파수( $f_{p e a k}$ )로 설정하고 감쇠 계수를 적용하여 Fig. 1과 같이 송신파형을 생성하였다(Algorithm 1). 이후 송신파형은 주파수 영역에서 적용되므로 고속 푸리에 변환(Fast Fourier Transform, FFT)을 통해 주파수 영역으로 변환하였다.

https://cdn.apub.kr/journalsite/sites/ggekseg/2026-029-02/N0540290205/images/ggekseg_2026_292_131_F1.jpg

Fig. 1

Generated time-domain Ricker wavelet.

Algorithm 1.

Damped Ricker Wavelet Generation.

Input:	$f_{m a x}$ , $t_{m a x}$ , Δt, γ
1.	Set $f_{p e a k}$ = $f_{m a x}$ / 3
2.	Define time axis t = 0 : Δt : $t_{m a x}$
3.	Set source delay time t₀ = 1 / $f_{p e a k}$
4.	For each time sample tᵢ do: source(tᵢ) = [1 − 2{π $f_{p e a k}$ (tᵢ − t₀)}²] exp[−{π $f_{p e a k}$ (tᵢ − t₀)}²]
5.	Apply temporal damping: source(tᵢ) ← source(tᵢ) exp(−γ tᵢ)
6.	Transform source wavelet to frequency domain: SRC(fⱼ) = RFFT{ $s o u r c e_{d a m p e d}$ (tᵢ) }

해석해를 계산하기 위해 식 (30)과 (31)의 피적분 함수를 모든 주파수, 파수 성분에 대해 계산한다. 그 결과로 $n_{k} \times n_{f}$ 크기를 갖는 두 2차원 배열이 계산되며 이를 각각 $u$ , $w$ 라 명명하였다.

영역 변환을 수행하기 전, 식 (28)과 (29)에서 언급된 것처럼 $e^{i k_{x} x_{0}}$ 를 곱하여 $x = x_{0}$ 에 위치하는 점하중이 작용하는 경우를 구현할 수 있다. 또한 주파수–파수 영역 해석해의 역 푸리에 변환 전에 시간영역 송신파형을 푸리에 변환하여 곱함으로써 통상적인 탄성파의 파형 특성(band-limited wavelet)을 갖는 파동장을 구할 수 있다. 임펄스 형태의 송신원을 디락–델타 함수로 표현하고 풀이한 해석해는 파동장의 그린함수(Green’s function)를 나타내기 때문에, 송신파형을 주파수 영역에서 그린함수와 곱해줌으로써 송신원 에너지 파형이 매질과 커플링되는 현상을 모사하게 된다(Algorithm 2).

주파수–파수 영역 파동장 $u$ 와 $w$ 를 이중 푸리에 변환하여 시간–공간 영역 파동장을 얻을 수 있다. 식 (30)과 (31)에서 관찰되는 바와 같이, 푸리에 변환의 부호 규약(sign convention)에 따라 시간축에서는 $- i ω t$ 를 정방향(forward)으로, $+ i ω t$ 를 역방향(inverse)으로 정의하며, 공간–파수 변환에서는 파동의 진행 방향을 양의 축 방향으로 정의하기 위해 부호를 달리 적용한다. 이산화된 자료의 빠른 푸리에 변환을 위해 고속 푸리에 변환(FFT) 알고리즘을 사용한다. Algorithm 3과 4는 앞서 계산된 피적분함수에 푸리에 변환을 적용하는 과정을 보여준다.

Algorithm 2.

Integrand Functions for Lamb’s Problem.

Input:	$n_{k}$ , $n_{x}$ , $n_{f}$ , α, β, μ, x₀, SRC
1.	Define wavenumber array: kᵢ = fftfreq( $n_{x}$ , Δx) (i = 1, ..., nₓ)
2.	For each frequency sample fⱼ (j = 1, ..., $n_{f}$ ): ωⱼ = 2π fⱼ Ωⱼ = ωⱼ − iγ
3.	Compute vertical decay terms (Re ≥ 0): p(kᵢ, fⱼ) = sqrt( kᵢ² − (Ωⱼ/α)² ) q(kᵢ, fⱼ) = sqrt( kᵢ² − (Ωⱼ/β)² )
4.	Compute Rayleigh denominator: A = 2kᵢ² − (Ωⱼ/β)² R = A² − 4kᵢ² p(kᵢ, fⱼ) q(kᵢ, fⱼ)
5.	Compute integrand functions with source position x₀: u(kᵢ, fⱼ) = −(i/μ) · kᵢ{A − 2p(kᵢ,fⱼ)q(kᵢ,fⱼ)} / R · exp(ikᵢx₀) w(kᵢ, fⱼ) = (1/μ) · p(kᵢ,fⱼ)(Ωⱼ/β)² / R · exp(ikᵢx₀)
6.	Remove DC component (k = 0): u(\|kᵢ\| < ε, :) = 0, w(\|kᵢ\| < ε, :) = 0
7.	Apply source wavelet spectrum: u(kᵢ, fⱼ) ← u(kᵢ, fⱼ) SRC(fⱼ) w(kᵢ, fⱼ) ← w(kᵢ, fⱼ) SRC(fⱼ)

Algorithm 3.

Spatial Fourier Transform and Cropping (k→x).

Input:	u, w, nₓ, x, $x_{c r o p}$
1.	Perform FFT along wavenumber axis (k→x): $u_{x f}$ = FFT{u} (axis = 0) $w_{x f}$ = FFT{w} (axis = 0)
2.	Use a sufficiently wide spatial model to suppress wrap-around artifacts: x = 0, Δx, ..., (nₓ − 1)Δx
3.	Crop the target receiver range after the spatial transform: ${u_{x f}}^{c r o p}$ = $u_{x f}$ [x ∈ $x_{c r o p}$ , :] ${w_{x f}}^{c r o p}$ = $w_{x f}$ [x ∈ $x_{c r o p}$ , :]

Algorithm 4.

Temporal Inverse Fourier Transform (f→t).

Input:	${u_{x f}}^{c r o p}$ , ${w_{x f}}^{c r o p}$ p, nₜ, t, γ
1.	Perform IRFFT along frequency axis (f→t): $u_{x t}$ = IRFFT{ ${u_{x f}}^{c r o p}$ } (axis = 1, n = nₜ) $w_{x t}$ = IRFFT{ ${w_{x f}}^{c r o p}$ } (axis = 1, n = nₜ)
2.	Take real part: $u_{x t}$ ← Re( $u_{x t}$ ), $w_{x t}$ ← Re( $w_{x t}$ )
3.	Apply damping compensation: $u_{x t}$ (x, t) ← $u_{x t}$ (x, t) exp(γt) $w_{x t}$ (x, t) ← $w_{x t}$ (x, t) exp(γt)
4.	Transpose to (t × x): $u_{p l o t}$ = $u_{x t}^{T}$ $w_{p l o t}$ = $w_{x t}^{T}$

지금까지 설명한 과정을 통해 2차원 반무한 매질에서 지표면에 수직 방향으로 점하중이 작용할 때 지표면에서 측정되는 변위를 계산하였다. 계산 결과인 변위 파동장 $u$ 와 $w$ 는 각각 Fig. 2(a), 2(b)와 같다. 매질의 P파 속도는 1000 m/s, S파 속도는 500 m/s, 밀도를 2000 kg/m³으로 설정하고, 최대 주파수 50 Hz 리커 함수를 송신 파형으로 사용하여 기록시간 1 초의 파동장을 계산하였다. 모델링 공간 격자 크기는 1 m이며, x축 최대 거리를 1500 m로 설정하였다. Fig. 2(a), 2(b)는 송신원의 x 축 위치를 750 m, Fig. 2(c), 2(d)는 200 m로 계산한 결과이다.

https://cdn.apub.kr/journalsite/sites/ggekseg/2026-029-02/N0540290205/images/ggekseg_2026_292_131_F2.jpg

Fig. 2

Wavefields calculated from the analytical solution of Lamb’s problem. (a), (b) Cases in which the source is located at the center of the model ( $x_{0} = 750 m$ ); (c), (d) cases in which the source is located near the model boundary ( $x_{0} = 200 m$ ).

Fig. 2는 탄성 반무한 공간에서 수직 방향의 점하중이 작용할 때 발생하는 파동장의 수직 및 수평 변위를 시각적으로 나타낸 것이다. 이러한 결과는 해석해 기반 계산 결과가 탄성파의 주요 전파 특성을 잘 재현함을 확인할 수 있다. Fig. 2(c), 2(d)에서는 FFT의 주기적 가정으로 인해 공간 방향 두루마리 현상이 발생함을 확인할 수 있다. 시간 방향 두루마리 효과는 복소 각주파수를 사용하여 감쇠시키는 것이 가능한 반면, 공간 방향으로 같은 방식의 감쇠를 적용하여 두루마리 현상을 완화하는 것은 순환 경계조건(cyclic boundary condition)의 구조 때문에 구현하기가 어렵다. 따라서 본 연구에서는 공간 영역을 관심 영역보다 충분히 넓게 설정한 뒤, 관심 영역만 절단(cropping)하는 방식을 사용하였다. 따라서 매질의 x방향을 관심 영역보다 넓은 4,096 m로 설정한 후, 계산된 파동장을 잘라내는(crop) 과정을 수행하였다(Fig. 3).

https://cdn.apub.kr/journalsite/sites/ggekseg/2026-029-02/N0540290205/images/ggekseg_2026_292_131_F3.jpg

Fig. 3

Suppression of spatial wrap-around artifacts by computing the wavefield over an extended computational domain and cropping only the target region of interest.

주파수–파수 영역 해석해 기반 탄성 파동장 계산 결과의 타당성을 검증하기 위해, 동일한 탄성 매질 조건과 송신원 위치를 사용하여 2차원 시간영역 탄성 유한요소법(Finite Element Method, FEM) 수치 모델링을 수행하였다. 수치 모델링에서는 자유면(free-surface) 경계조건을 적용한 2차원 등방성(isotropic) 탄성파동 방정식을 사용하였으며, 해석해 풀이와 동일하게 P파 속도는 1000 m/s, S파 속도는 500 m/s, 밀도를 2000 kg/m³으로 설정하였다. 또한 자유표면의 $x_{0} = 200 m$ 위치에 수직 방향의 점하중을 적용하여 해석해와 동일한 조건을 구현하였다.

Fig. 4는 Fig. 3의 해석해 결과와 FEM 기반 탄성파 수치 모델링 결과를 비교한 것이다. Fig. 4(a), 4(c)는 각각 FEM 기반 수평 및 수직 변위 결과를 나타내며, Fig. 4(b), 4(d)는 Fig. 3의 해석해 결과와 FEM 결과 사이의 잔차(residual)를 나타낸다. 잔차 단면에서 대부분의 영역이 매우 작은 진폭으로 나타나며, 주요 파동 이벤트의 도달시간과 파형 형상이 전반적으로 잘 일치하는 것을 확인할 수 있다. 이는 푸리에 영역 기반 Lamb 문제의 해석해 계산 과정이 자유표면을 포함하는 탄성파의 주요 전파 특성을 적절히 재현함을 의미한다.

https://cdn.apub.kr/journalsite/sites/ggekseg/2026-029-02/N0540290205/images/ggekseg_2026_292_131_F4.jpg

Fig. 4

Comparison between the analytical solution shown in Fig. 6 and elastic-wave numerical modeling results obtained using the finite-element method (FEM). (a) FEM horizontal displacement, (b) residual between Fig. 6(a) and the FEM horizontal displacement, (c) FEM vertical displacement, and (d) residual between Fig. 6(b) and the FEM vertical displacement.

몇 개의 수신기 위치에서 파동장을 추출하여 트레이스를 비교한 결과(Fig. 5), 해석해와 FEM 기반 수치 모델링 결과는 전체적인 도달시간과 주요 파형 특성에서 매우 잘 일치하는 것으로 나타났다. 특히 직접파와 표면파 이벤트의 도달 시점 및 위상 특성이 전반적으로 일치하였으며, 그 차이 역시 대부분의 구간에서 매우 작은 진폭으로 확인되었다. 또한 수평 및 수직 변위 성분 모두에서 해석해와 FEM 결과가 안정적으로 일치하는 것을 확인할 수 있다. 전체 파동장에 대한 상관계수를 계산한 결과, 수평 변위와 수직 변위 성분의 상관계수는 각각 0.9948, 0.9948로 나타났으며, 이를 통해 해석해와 FEM 기반 수치 모델링 결과가 매우 높은 수준으로 일치함을 정량적으로 확인할 수 있었다.

https://cdn.apub.kr/journalsite/sites/ggekseg/2026-029-02/N0540290205/images/ggekseg_2026_292_131_F5.jpg

Fig. 5

Trace comparison between the analytical solution and FEM-based numerical modeling results at multiple receiver locations. (a), (c), and (e) Horizontal displacement components; (b), (d), and (f) vertical displacement components.

결 론

이 논문에서는 Lamb 문제의 가장 단순한 경우를 가정하여, 탄성 반무한 공간에서 수직 점 하중에 의해 발생하는 주파수–파수 영역 파동장의 적분 형태 해를 유도하고 이를 이산 고속 푸리에 변환을 이용해 구현하는 과정을 설명하였다. 구현 과정에서 유의해야 할 주요 사항을 정리하면 다음과 같다.

1. 송신원 위치가 원점이 아닌 경우, 지표 응력 경계조건에서 중심이 원점이 아닌 디락–델타 함수를 사용한다. 이 경우 푸리에 영역에서 이동 정리(shifting theorem)를 적용하여 송신원 위치 변화는 복소 지수함수 형태의 위상항을 곱하는 방식으로 구현될 수 있다.

2. 역 푸리에 변환과정에서 시간 및 공간 방향의 두루마리 현상이 발생할 수 있는데, 시간 방향 두루마리 현상은 복소 각주파수를 이용해 감쇠시킬 수 있으나 공간 방향의 두루마리 현상은 순환 경계조건 때문에 동일한 방식을 적용하기 어렵다. 따라서 계산 영역을 관심 영역보다 충분히 넓게 설정한 뒤, 불필요한 파동장을 잘라내는 방식을 사용하는 것이 좋다.

3. Lamb 문제의 주파수–파수 영역 파동장에서 수평성분과 수직성분은 공간파수 축을 기준으로 각각 기함수와 우함수로 주어지므로 역 푸리에 변환을 위한 음수 축의 파동장을 대칭성으로 계산할 때 부호를 유의해야 한다.

이 논문에서는 해석적인 복소 적분 기법을 적용하여 닫힌 형태의 해를 구하는 복잡한 풀이 기법은 다루지 않았지만, 푸리에 변환 기반 구현 과정을 직관적으로 설명함으로써 기초적인 파동방정식 해를 이해할 수 있도록 하였다. 이를 바탕으로 이론 지진학과 탄성파 탐사를 처음 접하는 입문자들에게 탄성파 해석해와 파동장 모델링의 기본 개념을 이해하는 데 하나의 길잡이가 되기를 기대한다.

Acknowledgements

이 성과는 정부(과학기술정보통신부)의 재원으로 한국연구재단의 지원을 받아 수행된 연구임(RS-2025-24523155).

References

Aki, K. and Richards, P. G., 2002, Quantitative Seismology, 2nd ed., University Science Books.

Bard, P. Y. and Bouchon, M., 1985, The two-dimensional resonance of sediment-filled valleys, Bull. Seismol. Soc. Am., 75(2), 519-541.

10.1785/BSSA0750020519

Cagniard, L., 1962, Reflection and Refraction of Progressive Seismic Waves, McGraw-Hill.

Červený, V., 2001, Seismic Ray Theory, Cambridge University Press.

Chao, C. C., 1960, Dynamical response of an elastic half-space to tangential surface loading, J. Appl. Mech., 27(3), 559-567.

10.1115/1.3644041

De Hoop, A. T., 1960, A modification of Cagniard’s method for solving seismic pulse problems, Appl. Sci. Res., 8, 349-356.

10.1007/BF02920068

Emami, M. and Eskandari-Ghadi, M., 2019, Lamb’s problem: A brief history, Mathematics and Mechanics of Solids, 25(3), 501-514.

10.1177/1081286519883674

Feng, X. and Zhang, H., 2018, Exact closed-form solutions for Lamb’s problem, Geophys. J. Int., 214(1), 444-459.

10.1093/gji/ggy131

Feng, X. and Zhang, H., 2020a, Exact closed-form solutions for Lamb’s problem–II: A moving point load, Geophys. J. Int., 223(2), 1446-1459.

10.1093/gji/ggaa380

Feng, X. and Zhang, H., 2020b, Exact closed-form solutions for Lamb’s problem–III: The case for buried source and receiver, Geophys. J. Int., 224(1), 517-532.

10.1093/gji/ggaa485

Graff, K. F., 1975, Wave Motion in Elastic Solids, Dover Publications.

Hadley, D. M. and Helmberger, D. V., 1977, Generalized ray theory for a layered medium, Geophys. J. R. Astron. Soc., 51(3), 401-416.

Hu, S. and Jiang, X., 2024, A note on exact closed-form solutions for Lamb’s problem with an arbitrary Poisson’s ratio, Bull. Seismol. Soc. Am., 114(5), 2277-2286.

10.1785/0120240028

Johnson, L. R., 1974, Green’s function for Lamb’s problem, Geophys. J. R. Astron. Soc., 37, 99-131.

Kausel, E., 2013, Lamb’s problem at its simplest, Proc. R. Soc. A, 469(2149), 20120462.

10.1098/rspa.2012.0462

Lamb, H., 1904, On the propagation of tremors over the surface of an elastic solid, Phil. Trans. R. Soc. Lond. A, 203, 1-42.

10.1098/rsta.1904.0013

Lapwood, E. R., 1949, The disturbance due to a line source in a semi-infinite elastic medium, Phil. Trans. R. Soc. Lond. A, 242(841), 63-100.

10.1098/rsta.1949.0005

Miller, G. F. and Pursey, H., 1954, The field and radiation impedance of mechanical radiators on the free surface of a semi-infinite isotropic solid, Proceedings of the Royal Society of London A, 223(1155), 521-541.

10.1098/rspa.1954.0134

Mooney, H. M., 1974, Some numerical solutions for Lamb’s problem, Bull. Seismol. Soc. Am., 64(2), 473-491.

10.1785/BSSA0640020473

Pekeris, C. L., 1955a, The seismic surface pulse, Proc. Natl. Acad. Sci. U.S.A., 41, 469-480.

10.1073/pnas.41.7.46916589699PMC528118

Pekeris, C. L., 1955b, The seismic buried pulse, Proc. Natl. Acad. Sci. U.S.A., 41, 629-639.

10.1073/pnas.41.9.62916589718PMC528151

Pinney, E., 1954, Surface motion due to a point source in a semi-infinite elastic medium, Bull. Seismol. Soc. Am., 44(4), 571-596.

10.1785/BSSA0440040571

Rayleigh, L., 1885, On waves propagated along the plane surface of an elastic solid, Proceedings of the London Mathematical Society, 1, 4-11.

10.1112/plms/s1-17.1.4

Richards, P. G., 1979, Elementary solutions to Lamb’s problem for a point source and their relevance to three-dimensional studies of spontaneous crack propagation, Bull. Seismol. Soc. Am., 69(4), 947-956.

10.1785/BSSA0690040947

Geophysics and Geophysical Exploration ISSN:1229-1064(Print) 2384-051X(Online) 지구물리와 물리탐사

Preview

A Tutorial on the Analytical Solution of Lamb’s Problem and Its Numerical Implementation

ABSTRACT

MAIN

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

(24)

(25)

(26)

(27)

(28)

(29)

(30)

(31)

Table 1

Definitions and corresponding equations for the principal parameters used in the analytical solution calculation.

Fig. 1

Generated time-domain Ricker wavelet.

Algorithm 1.

Damped Ricker Wavelet Generation.

Algorithm 2.

Integrand Functions for Lamb’s Problem.

Algorithm 3.

Spatial Fourier Transform and Cropping (k→x).

Algorithm 4.

Temporal Inverse Fourier Transform (f→t).

Fig. 2

Wavefields calculated from the analytical solution of Lamb’s problem. (a), (b) Cases in which the source is located at the center of the model (x0=750m); (c), (d) cases in which the source is located near the model boundary (x0=200m).

Fig. 3

Suppression of spatial wrap-around artifacts by computing the wavefield over an extended computational domain and cropping only the target region of interest.

Fig. 4

Fig. 5

Trace comparison between the analytical solution and FEM-based numerical modeling results at multiple receiver locations. (a), (c), and (e) Horizontal displacement components; (b), (d), and (f) vertical displacement components.

Acknowledgements

References

Wavefields calculated from the analytical solution of Lamb’s problem. (a), (b) Cases in which the source is located at the center of the model ( $x_{0} = 750 m$ ); (c), (d) cases in which the source is located near the model boundary ( $x_{0} = 200 m$ ).