2020년 12월 6일 일요일

아벨의 적분 방정식(Abel's Integral Equation)

[경고] 아래 글을 읽지 않고 "아벨의 적분 방정식"을 보면 바보로 느껴질 수 있습니다.


대부분의 적분 방정식(integral equation)은 정확한 해법이 존재하지 않기 때문에 피적분 함수를 적절한 기저 함수(basis function)로 분해해서 수치 해석적으로 적분 방정식의 해를 근사화한다. 예외적으로 적분 방정식이 깔끔하게 풀리는 경우도 어쩌다 한 번 정도는 있다. 그 대표적인 예가 다음과 같은 아벨의 적분 방정식(Abel's integral equation)[1]이다.

                  (1)

여기서 $f(0)$ = $0$, $0 < \alpha < 1$이다. 아벨의 적분 방정식은 $t$ = $x$에서 피적분 함수가 발산하는 매우 특이한 성질을 가진 특이 적분 방정식(singular integral equation)이다. 식 (1)의 풀이법은 아벨Niels Henrik Abel(1802–1829)이 1823년아벨 21세, 조선 순조 시절에 제안했다. 베타 함수(beta function)를 적용하면 식 (1)을 쉽게 증명할 수 있다. 먼저 식 (1)의 왼쪽 식을 $(x-t)^{1- \alpha}$로 나눈 후 $t$에 대해 적분한다.

                  (2)

식 (2)에 나타난 무리 함수의 적분은 베타 함수로 표현되므로, 적분 결과는 $\tau$에 대해서 항상 상수가 된다.

                      (3)

따라서 식 (3)을 식 (2)에 대입한 후 $x$에 대해 미분하면 식 (1)이 증명된다.

                      (4)

식 (1)의 오른쪽 식에 나타난 적분에 미분 연산 $d/dx$를 적용하기 위해 부분 적분(integration by parts)을 기반으로 피적분 함수의 특이성을 제거한다.

                      (5)

여기서 $f'(t)$ = $df(t)/dt$이다. 식 (5) 우변의 첫째 항을 $x$에 대해 미분하면 다음 결과를 얻는다.

                      (6)

다음 단계로 둘째 항도 미분을 해야 하지만, 피적분 함수와 적분 구간에 $x$가 모두 있어서 미분이 어렵다. 이런 어려움을 피하기 위해 $(x-a)^\alpha$를 뉴턴의 이항 정리(Newton's binomial theorem)로 전개해서 적분하고 미분한다.

                      (7)

                      (8)

여기서 $(\alpha)_k$는 포흐하머 기호(Pochhammer symbol) 혹은 하강 계승(falling factorial)이다. 식 (8)을 더 쉽게 유도하려면 정적분의 미분 관계를 이용해야 한다. 식 (6)과 (8)을 합쳐서 $\phi(x)$를 더 간단하게 표현한다. 

                      (9)

식 (1)의 적분 구간을 다음처럼 $0$에서 무한대로 바꾸어서 조금 다른 공식을 만들 수도 있다.

             (10)

여기서 $\lim_{x \to \infty} f(x)$ = $0$이다. 식 (10)의 유도는 식 (2)와 비슷한 과정을 거친다.

                      (11)

식 (5)처럼 피적분 함수의 특이성을 없애기 위해 식 (10)의 오른쪽 식을 부분 적분한다.

                      (12)

다음 단계로 식 (12)의 우변을 $x$에 대해 미분하여 정리한다.

                      (13)

                      (14)

식 (13)과 (14)를 합친 후, 식 (10)의 오른쪽 식에 대입해서 최종 결과를 얻는다.

                      (15)

아벨의 적분 방정식에 나오는 $\phi(x)$는 어떤 함수든지 가능하므로, 식 (1)의 적분 방정식을 적분 변환(integral transform) 형태로 정의하기도 한다. 즉 아벨 변환(Abel transform)아벨 역변환(inverse Abel transform)이란 이름을 붙여서 식 (1)을 다음처럼 표기할 수 있다.

                      (16)

여기서 $df(\chi)$ = $f'(\chi) d\chi$이다. 식 (16)의 마지막식처럼 표현한 적분은 리만–스틸체스 적분(Riemann–Stieltjes integral)이라 부른다.
아벨의 적분 방정식은 분수 미적분학(fractional calculus)으로 개념화하기도 한다. 아벨 사후 3년 뒤인 1832년리우빌 23세, 조선 순조 시절에 리우빌Joseph Liouville(1809–1882)이 분수 미분과 적분(fractional differentiation and integration)을 제안했고, 1847년리만 21세, 조선 헌종 시절에는 리만Bernhard Riemann(1826–1866)이 분수 미적분학을 복소 영역(complex domain)으로 확장했다. 한참 시간이 흐른 1893년헤비사이드 43세, 조선 고종 시절에 헤비사이드Oliver Heaviside(1850–1925)연산 미적분학(operational calculus)을 새롭게 제안하면서 분수 미적분학을 다시 강조했다. 그래서 식 (1)에 나오는 적분을 다음과 같은 리만–리우빌 적분(Riemann–Liouville integral)이라고 부른다.

                  (17)

식 (17)의 정의를 이용해서 식 (1)의 왼쪽 식을 다시 쓴다.

                  (18)

순전히 연산자 관점으로만 식 (18)의 좌변을 본다. 그러면 분수 적분인 식 (18)을 다시 분수 미분해주면 원래 함수 $\phi(x)$가 당연히 나온다.

                  (19)

따라서 분수 미적분학을 이용해도 식 (1)과 같은 관계가 유도된다.


   1. 기본(basics)   

아벨의 적분 방정식은 아래처럼 다양한 방식으로 표현될 수 있다.

[기본 표현식]

             (1.1)

             (1.2)

여기서 $f(0)$ = $0$, $0 < \alpha < 1$이다.

[증명]
식 (1.1)에 있는 변수를 $u = t^2$과 $v = x^2$로 치환해서 식 (1)과 비교하여 증명한다. 비슷하게 식 (1.2)는 변수 치환을 통해 식 (10)으로 만든다.
______________________________

라돈 변환(Radon transform)처럼 원통 좌표계에 등장하는 적분 방정식의 해를 구할 때에는 식 (1.1), (1.2)가 매우 편리하다.


[참고문헌]
[1] B. V. Khvedelidze, "Abel integral equation," Encyclopedia of Mathematics, 2011.

[다음 읽을거리]

2020년 11월 28일 토요일

라돈 변환(Radon Transform)

[경고] 아래 글을 읽지 않고 "라돈 변환"을 보면 바보로 느껴질 수 있습니다.


[그림 1] 라돈 변환과 컴퓨터 단층 촬영(출처: wikipedia.org)

푸리에 급수(Fourier series)푸리에 변환(Fourier transform)현대 수학을 만든 시초라 할 수 있다. 푸리에 변환에 버금가는 대단한 기여를 한 적분 변환중 하나는 라돈 변환(Radon transform)이다. 라돈 변환은 수학자 라돈Johann Radon(1887–1956)이 1917년라돈 30세, 일제 식민지 시절에 제안했다[1]. 사영(projection)을 이용해 원래 물체의 영상을 복원하는 라돈 변환의 대표적 응용은 [그림 1]에 소개한 컴퓨터 단층 촬영(computed tomography, CT)이다. CT 기술은 EMI의 공학자인 하운스필드Godfrey Hounsfield(1919–2004)에 의해 1972년하운스필드 53세, 박정희 정부 시절에 최초로 상용화되었다.

[그림 2] 미국을 방문한 비틀즈(출처: wikipedia.org)

음반 회사로 유명한 EMI는 비틀즈(The Beatles) 때문에 엄청난 매출과 수입을 올렸다. 비틀즈가 번 돈은 EMI가 만들던 CT 기술에도 투입되어서 인류를 위한 의료 기술 개발에 훌륭하게 쓰였다. 비틀즈의 노래를 사랑한 팬들은 알게 모르게 CT 기술에 대한 간접 지원을 한 것이었다. 돈을 매개로 한 예술과 과학의 놀라운 만남이다.

[그림 3] CT 주사기(scanner)의 내부 구조(출처: wikipedia.org)

[그림 3]은 전형적인 CT 주사기(scanner)의 내부 구조를 보여준다. X선관(X-ray tube) T에서 나온 X선 X는 공동[빈 구멍] 위치에 놓일 물체에 일부 흡수되면서 X선 검출기(X-ray detector) D까지 거의 직선으로 진행한다. 더 세부적으로 보면, X선은 직진성이 매우 높아서 송신기에서 수신기까지 가는 경로는 거의 직선으로 가정할 수 있고, 물체에 의해 흡수되기 때문에 물체의 종류에 따라 투과율이 달라진다. 예를 들어 전기장의 거리별 감쇠(attenuation)는 $|E(z)|$ = $|E_0| e^{-\alpha z}$로 표현되므로, [그림 3, 4]처럼 X선관에서 발생한 X선이 물체를 거쳐 검출기에 수신되는 복사 선속(radiant flux) 혹은 복사 전력(radiant power) $I_\theta(s)$는 다음처럼 공식화한다.

                  (1)

여기서 $I_0$는 X선관의 복사 선속, $\alpha(s, t; \theta)$는 물체를 통과할 때 생기는 위치별 감쇠 상수, $p_\theta(s)$는 사영 함수(projection function)이다. 식 (1)을 정리해서 사영 함수를 물체의 감쇠 상수로 표기할 수도 있다.

                  (2)

여기서 $I_\theta(s)/I_0$는 X선 검출기의 수신 전력 비율이다. 따라서 라돈 변환의 핵심인 사영 함수는 복사 선속에 대한 로그 함수 비율과 등가이다. 즉 라돈이 제안한 라돈 변환은 [그림 3]에 보여준 물체에 대한 X선(X-ray) 투과와 밀접히 연결된다.

[그림 4] 라돈 변환의 좌표계(출처: wikipedia.org)
(1) 물체(object)
(2) 광선의 시작선(starting line of rays)
(3) 광선의 종료선(ending line of rays)
(4) 투과 광선(transmission ray)
(5) 자료 원(datum circle)
(6) 원점(origin)
(7) 사영 함수(projection function) $p_\theta(s)$

[그림 4]는 라돈 변환을 정의할 때 사용하는 좌표계를 표현한다. [그림 3, 4]를 비교하면 수학적인 라돈 변환을 물리적 실재로 구현한 결과가 CT 기술임을 잘 알 수 있다. [그림 4]에서 주어진[혹은 고정된] $s, \theta$에 대해 점 $(x, y)$를 새로운 좌표축 $t$로 표현하면 다음과 같다.

                  (3)

여기서 $s$ = $x \cos \theta + y \sin \theta$, $t$ = $-x \sin \theta + y \cos \theta$, $x^2 + y^2$ = $s^2 + t^2$이다. 식 (3)을 간편하게 2차원 회전 행렬(rotation matrix)처럼 쓸 수도 있다.

                  (4)

여기서 점 $(s, t)$를 위한 좌표축 $s, t$는 기준 좌표축 $x, y$를 $\theta$만큼 회전시켜 얻는다. 시작선에서 나와서 종료선까지 적분하는 사영 함수(projection function) $P_\theta(s)$는 다음과 같은 선 적분(line integral)으로 정의한다.

                  (5)

여기서 $f(x, y)$는 고정점 $z$ = $z_0$에 대한 물체의 2차원 절편(slice)이며, $(x, y)$는 식 (3)의 정의를 사용한다. 사영 함수 $p_\theta(s)$는 공동도(空洞圖) 혹은 사이노그램(sinogram)이라고도 한다. 공동도는 구멍(sinus, opening)에 대한 그림(gram)이란 뜻이다. 변수 $t$가 변하는 범위를 전체 실수로 확장하면 식 (5)는 라돈 변환이 된다.

                  (6)

각도를 $\theta + \pi$로 바꾸면, $p_{\theta + \pi}(s)$ = $p_\theta(-s)$ 혹은 $p_\theta(s)$ = $p_{\theta + \pi}(-s)$가 항상 성립한다.

[그림 5] 푸리에 절편 정리(출처: wikipedia.org)

사영 함수는 푸리에 변환 관점으로 관찰해야 숨겨진 내면을 볼 수 있다. 푸리에 절편 정리(Fourier slice theorem) 혹은 사영-절편 정리(projection-slice theorem)에 따라 사영 함수의 푸리에 변환을 구한다. 2차원에 대한 푸리에 절편 정리는 다음과 같다. 

                  (7)

식 (7)에 제시한 푸리에 절편 정리에 의해, 사영 함수 $p(x)$의 푸리에 변환 $\mathfrak{F}[p(x)]$는 차원 하나가 축소된 푸리에 변환 $P(\xi)$ = $F(\xi, 0)$을 만든다. 이 개념을 $N$차원 공간으로 일반화해서 식 (7)을 $N$차원 푸리에 절편 정리로 확장한다.

                  (8)

여기서 $p(\cdot)$는 좌표 성분 $x_{N-m+1}, \cdots, x_N$에 대한 공간 영역의 사영(projection), $\mathfrak{F}_{N-m}[\cdot]$는 $x_1, \cdots, x_{N-m}$에 대한 푸리에 변환, $F(\cdot)$는 $\xi_{N-m+1}, \cdots, \xi_N$ 성분을 $0$으로 바꾼 파수 영역의 절편(slice)이다.
식 (7)에 의해 공간 영역(spatial domain)에서 회전 각도를 $\theta$만큼 돌리면서 구한 사영 함수 $p_\theta(s)$의 푸리에 변환 $P_\theta (\kappa)$는 파수 영역(spectral domain)에서도 각도 $\theta$만큼 회전한다. 따라서 다음 관계식을 보간(interpolation)해서 $f(x, y)$의 2차원 푸리에 변환 $F(\xi, \eta)$를 근사한다.

                  (9)

2차원 푸리에 절편 정리인 식 (9)를 이용해서 파수 영역의 2차원 푸리에 변환 $F(\xi, \eta)$를 쉽게 구할 수 있다고 오해할 수 있다. 하지만 파수 영역의 절편 함수(slice function) $P_\theta (\kappa)$는 원점 $(\xi, \eta)$ = $(0, 0)$를 지나는 직선이라서, 원점 부근의 표본[낮은 파수 성분]이 과대 평가되고 원점을 벗어난 표본[높은 파수 성분]은 과소 평가되는 심각한 문제가 있다.

[그림 6] 복소 영역에서 절편 함수의 분포(출처: wikipedia.org)

파수 영역의 절편 함수[파란색 선]가 가진 한계를 [그림 6]이 분명히 보여준다. 원점에는 직선이 잘 모이지만, 반지름이 커지면 직선간의 간격이 넓어져서 문제가 된다. 그래서 [그림 6]과 같은 파수 영역의 절편 함수를 적절히 보간(interpolation)해서 근사적으로 절편 표본간의 간격을 비슷하게 만든다. 그 다음에 보간한 $F(\xi, \eta)$의 푸리에 역변환을 취해서 물체의 2차원 영상 $f(x, y)$를 다음처럼 근사적으로 복원한다.

                  (10)

여기서 $\operatorname{Int}[{\bf A}]$는 집합(set) $\bf A$의 원소로 만드는 보간을 의미한다.
식 (10)은 근사가 많이 들어간 영상 복원법이라서, 수학적으로 더 엄밀하고 세련되게 이론을 전개할 필요가 있다. 이를 위해 원통 좌표계에서 사용하는 푸리에 변환인 한켈 변환(Hankel transform)을 식 (10)에 적용해서 정확한 라돈 역변환(inverse Radon transform)을 정의한다. 모든 각도 $\theta$에 대해 식 (9)를 계산한 후 복원한 물체의 영상 $f(x, y)$는 다음과 같다[3].

                  (11)

사영 함수의 푸리에 변환 $P_\theta (\kappa)$는 $\theta$에 대해 주기적이므로 푸리에 급수로 전개한다. 이 결과를 식 (11)에 넣고 한켈 변환처럼 정리해서 라돈 역변환을 얻는다.

                  (12)

여기서 $x$ = $\rho \cos \phi$, $y$ = $\rho \sin \phi$, $J_n (\cdot)$는 제$n$차 제1종 베셀 함수(Bessel function of the first kind)이다. 파수 영역에서 정의된 $P_n(\kappa)$는 $\theta$에 대한 평균값으로 정확히 구한다.

                  (13)

한켈 변환으로 라돈 역변환을 유도할 수 있지만, 권위 있는 원류는 라돈이 제안한 방법[1]이다. 방사선 원소인 라돈(Radon, Rn) 혹은 라듐(radium, Ra)과 이름이 비슷해서 위험하다는 오해를 받을 때도 있지만, 라돈은 인류에 기여한 아름다운 수학자이다. 라돈이 제안한 라돈 역변환은 사영 함수 $p_\theta(s)$로부터 시작한다. 디랙 델타 함수(Dirac delta function)를 도입해서 사영 함수를 2차원 공간에서 다시 표현한다.

                  (14)

[그림 4]에서 보여준 좌표축의 회전 각도 $\theta$에 대한 평균 사영 함수 $\bar p(s)$는 다음과 같다.

                  (15)

데카르트 좌표계(Cartesian coordinate system)를 원통 좌표계(circular cylindrical coordinate system)로 바꾸어서 디랙 델타 함수의 적분을 구한다.

                  (16)

여기서 $x$ = $\rho \cos \phi$, $y$ = $\rho \sin \phi$, $\rho > |s|$, $\theta_1$과 $\theta_2$는 $\delta(\cdot)$의 입력 변수에 대한 근이며 $\rho  \cos (\theta_{1,2} - \phi)$ = $s$가 성립한다. 식 (16)을 식 (15)에 넣고 $f(x, y)$의 평균값 $\bar f(\rho)$가 피적분 함수에 나타나도록 정리한다.

                  (17)

식 (17)은 아벨의 적분 방정식(Abel's integral equation)이므로, 다음과 같은 해석적인 해법으로 쉽게 해결된다.

          (18)

여기서 $\alpha$ = $1/2$이다. 따라서 $\bar f(\rho)$는 다음처럼 공식화된다.

                  (19)

여기서 식 (18)의 오른쪽 식에 나온 극한은 $0$이라 가정한다. 결국 반지름 $\rho$를 $0$으로 보내면, 함수값 $f(0, 0)$를 사영 함수의 적분으로 정확히 표현할 수 있다.

                  (20)

[그림 7] 일반적인 라돈 역변환을 위한 좌표 변환

임의의 위치에서 $f(x, y)$를 완벽하게 복원하기 위해, [그림 7]처럼 원점 $(0, 0)$을 $(x_0, y_0)$으로 보내는 좌표 변환에 따라 평균값 $\bar f(\rho)$와 $\bar p(s)$를 다음과 같이 일반화한다.

                  (21)

여기서 $\bar f(0, 0; r)$ = $\bar f(r)$, $\bar p(0, 0; q)$ = $\bar p(q)$, $q$는 $(x_0, y_0)$를 기준으로 양수나 음수가 된다. 식 (15), (17)에 식 (21)을 적용해서 식 (19)를 변형한다.

                  (22)

최종적으로 반지름 $r$이 $0$으로 가는 극한에 의해 식 (22)는 물체의 2차원 영상 $f(x, y)$가 된다.

                  (23)

식 (23)의 적분에는 함수 미분소가 있어서 다소 복잡하므로, 부분 적분을 적용해서 다음처럼 단순화한다.

                  (24)

[라돈 역변환의 조건]
라돈 역변환이 성립하기 위한 조건은 다음과 같다.
  • 물체의 영상 $f(x, y)$는 연속이다.
  • 식 (17)의 성립을 위해 다음 이중 적분은 수렴한다: 

  • 식 (22)를 만족하기 위해 임의의 점 $(x, y)$에서 다음 극한이 성립한다: $\lim_{r \to \infty} \bar f(x, y; r)$ = $0$ 

[참고문헌]
[1] J. Radon, "Über die Bestimmung von Funktionen durch ihre Integralwerte längs gewisser Mannigfaltigkeiten (On the determination of functions by their integral values along certain manifolds)," Berichte über die Verhandlungen der Königlich-Sächsischen Akademie der Wissenschaften zu Leipzig (Reports on the Proceedings of the Royal Saxonian Academy of Sciences in Leipzig), vol. 69, pp. 262–277, Apr. 1917. (In German, 방문일 2020-11-28)
[2] 박미경, 이의택, "Computed Tomography에 의한 절편 영상의 추정", 전자통신동향분석, 제12권, 제6호, pp. 50–57, 1997년 12월.
[3] E. Won, "Derivation of the inverse Radon transformation," Eunil Won's Home Page, Nov. 2007. (방문일 2020-12-05)