2022년 8월 20일 토요일

전자기 광선 추적(Electromagnetic Ray Tracing)

[경고] 아래 글을 읽지 않고 "전자기 광선 추적"을 보면 바보로 느껴질 수 있습니다.

[확인] 본 페이지는 exp(-iωt) 시간 약속을 사용하고 있습니다.


[그림 1] 전자기 광선 추적을 위한 좌표계

[표 1] 전자기 광선 추적에 쓰는 광선의 자료 구조, PRD(per-ray data)
이름
(Name)
자료형
(Data type)
의미
(Meaning)
전기장 계수(electric-field coefficient)complex3광선이 가진 전기장 계수의 3차원 벡터는 $\bar E_g$ = $(E_{gx}, E_{gy}, E_{gz})$로 설정;기하 광학장(geometrical optics field)의 계수를 의미; 각 성분은 페이저(phasor)를 사용하여 경로 길이 $s$ = $0$인 경우의 전기장 크기와 위상을 표현
수신 전기장(received electric field)complex3전기장 계수와 경로 길이까지 모두 고려한 실제로 수신되는 전기장 벡터 $\bar E_\text{rx}$; 수신기에서 감지되는 모든 송신 광선을 포함한 값
경로 길이(path length)float 혹은 double광선이 이동한 총길이 $s$를 담는 변수; 기하 광학에 나오는 $A(s) e^{ik_0 s}$에 사용
반사 회수(# of reflection)unsigned int광선이 물체면에서 반사되는 회수를 저장; 반사와 회절 회수가 어느 이상이 되면 광선을 종료함
회절 회수(# of diffraction)unsigned int광선이 물체 모서리에서 회절되는 회수를 저장; 회절된 전기장은 약하기 때문에 회절 회수는 대략 $1$회만 생각함

전자파 신호의 송수신 특성을 계산하는 전자기 광선 추적(electromagnetic ray tracing)기존 광선 추적과 비슷하지만, 블린–퐁 반사 모형(Blinn–Phong reflection model)을 쓰지 않고 전자파의 반사(reflection)와 회절(diffraction)을 가능한 한 정확하게 적용한다. CG용 광선 추적을 채용하여 전자파의 산란 현상을 매우 빠르게 계산하는 전자기 광선 추적은 주로 통신 시스템의 전파 환경 분석에 사용된다[1]–[9]. 다만 [그림 1]처럼 전자기 광선 추적은 기존 광선 추적과는 다르게 송신기(transmitter, Tx)부터 출발해서 수신기(receiver, Rx)를 찾는다. 이러한 방식은 수신기부터 출발해 광선을 수동적으로 추적하는 광선 추적(ray tracing)은 아니고 송신기에서 광선을 능동적으로 던지는 광선 투사(光線投射, ray casting)가 맞는 용어이다. 하지만 광선 추적이 너무 유명해서 송신기에서 광선을 투사하는 방식도 어색하게나마 광선 추적이라 이름 붙인다. 또한 기존 광선 추적처럼 수신기에 들어오는 광선만 따라가면, 수신 전자파의 편파(polarization)를 고려하기 매우 어렵다. 왜냐하면 수신기로 들어오는 편파는 2종류가 있고 어느 편파가 수신될지 알 수 없어서 모든 편파를 다루어야 하기 때문이다. 즉, 송신기의 광선 투사보다 수신기의 광선 추적은 저장 공간과 정확한 산란 계산에 불리하다.
전자기 광선 추적의 공식화를 위해, 광선이 전파되는 특성은 기하 광학(geometrical optics, GO)의 관계식 $\bar E(s)$를 그대로 채택한다.

                  (1)

여기서 $s$는 전자파가 진행한 경로 길이(path length), 구면파인 경우 확산 인자(spreading factor)는 $A(s)$ = $\rho/(\rho + s)$, $\rho$는 구면파의 곡률 반경, $k$는 파수(wavenumber)이다. 이동 통신(mobile communication)에서 사용하는 안테나는 보통 수직 편파(vertical polarization)를 사용하므로, 전자기 광선 추적의 편파는 건물 벽에 전기장이 평행한 TE(횡전기 혹은 가로 전기, Transverse Electric)파를 주로 선택한다. 그래서 편파를 TE파로 고정하면 식 (1)의 벡터 $\bar E(s)$ 대신 스칼라 $E(s)$를 쓸 수도 있다. 임의의 안테나 편파를 다룰 때는 TE파와 TM(횡자기 혹은 가로 자기, Transverse Magnetic)파를 함께 사용한다. 프레넬 방정식(Fresnel equation)에 따라 TE파와 TM파가 임의 매질에서 반사되는 계수는 다음과 같다.

                  (2: TE파의 반사 계수)

                  (3: TM파의 반사 계수)

광선 추적에 필요한 반사 계수 $r_s$와 $r_p$는 입사 각도 $\theta_i$와 물체의 복소 유전 상수(complex dielectric constant) $\epsilon_r$로만 결정된다. 이상을 바탕으로 전자기 광선 추적을 OptiX 광선 추적 엔진(ray tracing engine)에서 실행하는 방법을 소개한다.

  • 광선 발생 프로그램(Ray Generation Program)
송신 안테나의 안테나 이득(antenna gain)과 관련된 벡터식 복사 패턴(vectorial radiation pattern)을 이용해 송신기에서 광선을 발생시킨다. 송신 안테나에서 복사되는 전기장 $\bar E(\bar r)$은 다음처럼 다양하게 공식화된다.

                  (4)

여기서 $V_\text{in}$은 안테나의 입력 전압, $\bar E_p(\theta, \phi)$ = $E_{p \theta}(\theta, \phi) \hat \theta + E_{p \phi}(\theta, \phi) \hat \phi$는 벡터식 복사 패턴(vectorial radiation pattern), $\bar E_p(\theta, \phi)$의 크기 제곱은 흔히 사용하는 안테나 이득(antenna gain) $G(\theta, \phi)$ = $|\bar E_p(\theta, \phi)|^2$이 된다. [표 1]에 제시한 PRD(per-ray data)가 가진 제$q$번 광선의 전기장 계수는 $\bar E_g$ = $\bar E_p(\theta_q, \phi_q)$로 둔다. 전기장 계수의 자료형인 complex3는 complex의 3차원 벡터 형태이다. 광선 발생 프로그램에서 광선은 물체를 고려하지 않고 모든 방향으로 골고루 혹은 무작위로 발사된다. 목적없이 임의로 발사하여 광선을 무작위로 추적하는 방식은 억지 광선 추적(brute-force ray tracing) 혹은 탈진된 광선 추적(exhausted ray tracing)이라 부른다. 여기서 제$q$번 광선이 발사되는 방향은 $(\theta_q, \phi_q)$로 생각하고, 발사하는 총광선 개수는 $Q$개라고 가정한다.

  • 최근접 명중 프로그램(Closest-hit Program)
광선이 물체와 최근접 명중(closest hit)하면, 이 물체에서 수신 안테나가 보이는지 확인한다. 수신 안테나가 보이면 TE파와 TM파의 반사 계수인 $r_s$와 $r_p$를 광선이 가진 전기장 계수 $\bar E_g$에 곱하여 반사파를 만들고, 송신 안테나에서 수신 안테나까지의 경로 길이 $s$를 대입해서 전기장 $\bar E(s)$를 근사적으로 결정한다. 쉬운 공식화를 위해 TE파와 TM파의 반사 계수를 다이애드(dyad)로 확장한 다이애드 반사 계수(dyadic reflection coefficient) $\bar{\bar{r}}_E$를 도입한다.

                  (5)

여기서 $\hat e_\parallel^i$와 $\hat e_\perp$는 각각 전기장 기준으로 설정한 입사 평면에 평행 및 수직인 단위 벡터, $\hat e_\parallel^r$은 전기장 기준 반사 평면에 평행한 단위 벡터이다. 그러면 송신 안테나에서 복사된 제$q$번 광선이 수신 안테나에 형성하는 전기장 $\bar E_q(s_q)$와 수신 전압(received voltage) $V_q(s_q)$는 다음과 같다.

                  (6a)

                  (6b)

여기서 $(\theta_q, \phi_q)$와 $(\vartheta_q, \varphi_q)$는 각각 송신과 수신 안테나에서 물체를 보는 제$q$번 광선의 각도, $s_q$는 제$q$번 광선이 송신 안테나에서 수신 안테나까지 가는 경로 길이, $\bar L_e(\theta, \phi)$는 수신 안테나의 벡터식 유효 길이(vectorial effective length), 수신 전압을 최대로 얻기 위해 $\bar L_e^*(\theta, \phi)$를 쓰는 켤레 정합(conjugate matching)을 한다. 건물 벽의 유전 상수와 반사 계수가 변하는 특성은 [6], [10]을 참고할 수 있다. 벡터식 유효 길이 $\bar L_e(\theta, \phi)$의 크기 제곱은 수신 안테나의 유효 면적(effective area) $A_e(\theta, \phi)$ = $|\bar L_e(\theta, \phi)|^2$이 된다. 또한 수신 안테나의 유효 길이 $\bar L_e(\theta, \phi)$는 송신 안테나의 벡터식 복사 패턴 $\bar E_p(\theta, \phi)$와 선형적으로 비례한다.

                  (7)

여기서 왼쪽과 오른쪽 식은 각각 전기장과 전력의 관계를 가진다.

[그림 2] 전자기 광선 추적에 쓰는 벡터의 정의

다만 식 (6)은 송신 광선(Tx ray)이 물체에서 반사되어 생기는 방향인 $\hat R$에서만 성립한다. [그림 2]처럼 수신기가 $\hat R$이 아닌 $\hat V$에 위치하면, 식 (5)보다 수신 전기장을 축소시켜야 한다. 그래서 퐁 반사 모형(Phong reflection model) 혹은 잘 알려진 포인팅의 정리(Poynting's theorem)에 따라, $\hat R$에서 벗어난 정도인 $\hat R \cdot \hat V$에 수신 전기장이 비례한다고 가정하여 식 (6)을 변형한다.

                  (8)

여기서 $\alpha$는 수신 비율을 근사화하는 물체의 반짝임 상수(shininess constant)이다. 반사 광선(reflection ray)의 방향 $\hat R$은 $\hat n, \hat T$로 딱 정해진다.

                  (9)

반짝임 상수 $\alpha$를 잘 모르면 단순하게 0.5로 둔다.[이렇게 두면 수신 전력이 내적 $\hat R \cdot \hat V$에 선형 비례한다.] 그 다음 단계로 물체에 의한 반사 광선을 다시 생성해서 발사하며, 이 반사 광선은 마치 송신 광선처럼 재귀적으로 처리한다. 이와 같이 제$q$번 반사 광선이 다시 반사 전기장 $\bar E_{rq}(s_{rq})$를 도출해서 우리가 얻기 원하는 수신 전기장 $\bar E_\text{rx}$에 모두 합산시킨다.

                  (10)

여기서 제$q$번 반사파의 경로 길이 $s_{rq}$는 송신 안테나에서 현재 물체까지 오는 경로 길이 $t_q$를 시작값으로 해서 반사가 생길 때마다 경로 길이를 계속 축적한다. 쉽게 말해서 식 (8)에 유도한 송신 광선의 수신 전기장 $\bar E_q(s_q)$처럼, [그림 1,2]에 제시한 반사 광선이 다시 물체와 만나서 수신 안테나를 볼 수 있으면, 수신 전기장을 $\bar E_{rq}(s_{rq})$만큼 만든다고 생각한다.


[그림 3] 반사 광선의 시작점을 표현하는 영상 전원

(a) 송신기에 의한 제1차 영상 및 회절 전원

(b) 제1차 전원에 의한 제2차 영상 및 회절 전원
[그림 4] 건물 벽에 의해 생성된 영상 및 회절 전원(출처: [6])

식 (8)은 기하 광학에 기반을 두어서 좋은 근사이기는 하지만, 퐁 반사 모형을 쓰기 때문에 정말 물리적으로 타당한지는 의문이다. 내적 $\hat R \cdot \hat V$를 없애고 물리 모형을 더 정확하게 표현하기 위해, [그림 3]에 제시한 영상 전원(image source)을 써서 전자기 광선 추적을 구현한다[3]–[9]. 영상 전원은 반사 광선이 시작되는 점으로 가정하므로, [그림 3]처럼 실제 전원의 위치 $A$를 경계면에 대칭해서 영상 전원을 위치 $B$에 정의한다. 영상 전원이 완벽히 만들어지기 위해서는 무한히 넓은 평면 구조가 필요하지만, 파장에 비해 건물 벽이 충분히 크다고 가정해서 건물 벽의 반대편에 영상 전원이 근사적으로 형성된다고 생각한다. 그러면 반사 광선이 $B$에서 나온다고 생각해서, 수신 안테나가 어디에 있든지 $\hat R \cdot \hat V$를 쓰지 않고 전기장을 계산할 수 있다. 대신 CG용 광선 추적과는 다르게 영상 전원을 중심으로 자료 구조를 구성해야 한다[3]–[6]. 예를 들어, 송신기(transmitter, Tx)가 만드는 전자파에 의한 건물 벽의 영상 전원을 보여주는 [그림 4(a)]를 고려한다. 송신기를 마주하는 벽은 2면이 있으므로, 영상 전원도 2개를 구성한다. 즉, 송신기와 마주보는 벽 #4에 의한 영상 전원은 #1[그림 4(a)에서 tube1으로 표기]이고, 벽 #5는 영상 전원 #2[그림 4(a)에서 tube2]를 등가적으로 발생시킨다.
영상 전원을 사용함으로써 얻어지는 또 하나의 장점은 회절(回折, diffraction) 처리에 있다. 손실 금속에 대한 회절 계수(diffraction coefficient)[11]까지 이미 근사적으로 공식화되었기 때문에, 전자파의 회절 계산은 간단한 연습 문제가 된다. 하지만 전자기 광선 추적에서 회절을 계산하기에는 엄청난 난관이 숨어있다. 전자파가 모서리(edge)에 입사하면 당연하게 회절이 생기지만, 모든 모서리마다 회절 계산용 광선이 들어와야 해서 광선 추적을 재귀 호출할수록 사용하는 광선이 기하 급수적으로 증가해서 계산 불능에 쉽게 빠진다. 이 문제를 영상 전원은 쉽게 해결할 수 있다. 영상 전원에 바탕을 둔 전자기 광선 추적은 억지 광선 추적(brute-force ray tracing)을 하지 않고, 건물 벽에 의해 형성된 영상 전원만 다루므로 회절을 고려해도 광선 개수가 기하 급수적으로 증가하지 않는다. 쉽게 말해 모서리가 있더라도 송신기 혹은 영상 전원에서 모서리로 가는 광선만 처리하므로, 회절까지 포함한 광선 추적의 계산량이 아주 크게 늘어나지 않는다. 이 개념을 이해하기 위해 [그림 4(a)]를 다시 본다. 송신기에서 모서리 a, d, e로 가는 광선은 회절을 만드는 원천인 회절 전원(diffraction source) #3, #4, #5[그림 4(a)에서 tube3, tube4, tube4로 표시]를 각각 생성한다. 이 회절 전원을 다시 시작점으로 삼아서 전자파의 전파 환경을 분석한다. 또한 [그림 4(b)]처럼 송신기가 만드는 제1차 전원은 다시 새로운 원천이 되어서 건물 벽과 모서리에 의한 제2차 영상 및 회절 전원을 재귀적으로 만든다. 예를 들어, 영상 전원 #2는 벽 #3에서 반사되므로 제2차 영상 전원 #6을 도출한다. 비슷하게 영상 전원 #1은 벽 #5에 따라 제2차 영상 전원 #7를 생성시킨다. 추가로 영상 전원 #2는 모서리 a, d에 제2차 회절 전원 #8, #9도 만든다.

(a) 건물 사이에 존재하는 송신기와 수신기

(b) 광선 추적을 위한 자료 구조 예시
[그림 5] 영상 전원을 쓰는 전자기 광선 추적의 자료 구조 구성 방식(출처: [3])

영상 및 회절 전원을 도입한 전자기 광선 추적을 위한 자료 구조(data structure)의 구성 방식을 이해하기 위해 [그림 5]를 꼼꼼히 관찰한다. 먼저 영상이나 회절 전원이 생기기 위해서는 광선 시작점에서 건물 벽이나 모서리가 보여야 한다. [그림 5(a)]의 송신기 관점에서 비출 수 있는 벽이나 모서리는 벽 #2, #16, 모서리 #1, #3, #9이다. 이 위치를 [그림 5(b)]와 같이 트리 혹은 나무(tree) 자료 구조의 마디(node)로 설정한다. 그 다음에 각 마디에 대해서 다시 재귀적으로 광선이 도달할 수 있는 벽이나 모서리를 추적한다. 모서리 #1는 모서리 #9와 벽 #16, 벽 #2도 모서리 #9와 벽 #16를 볼 수 있다. 모서리 #3은 더 많은 영역에 광선을 전파할 수 있어서 모서리 #9, #16, 벽 #16, 수신기(receiver, Rx)까지 전파를 전송한다. 이를 계속 확장해서 광선 추적용 자료 구조를 형성한 최종 결과는 [그림 5(b)]에 정리되어 있다. 물론 건물 벽이 반사하는 광선은 [그림 3, 4]와 같은 영상 전원으로 처리하고, 회절되는 광선은 회절 전원으로 간주하여 계산한다.

[참고문헌]
[1] K. J. Gladstone and J. P. McGeehan, "Computer simulation of multipath fading in the land mobile radio environment," IEE Proc. G - Electron. Circuits Syst., vol. 127, no. 6, p. 323–330, Dec. 1980.
[2] K. R. Schaubach, N. J. Davis and T. S. Rappaport, "A ray tracing method for predicting path loss and delay spread in microcellular environments," 42nd IEEE Veh. Technol. Soc. Conf.—Frontiers of Technology, vol. 2, pp. 932–935, May 1992.
[3] M. G. Sánchez, L. de Haro, A. G. Pino, and M. Calvo, "Exhaustive ray tracing algorithm for microcellular propagation prediction models," Electron. Lett., vol. 32, no. 7, pp. 624–625, Mar. 1996.
[4] M. F. Catedra, J. Perez, F. Saez de Adana, and O. Gutierrez, "Efficient ray-tracing techniques for three-dimensional analyses of propagation in mobile communications: application to picocell and microcell scenarios," IEEE Antennas Propag. Mag., vol. 40, no. 2, pp. 15–28, Apr. 1998.
[5] 손해원, 명로훈, "광선 추적법을 이용한 마이크로셀 전파환경 예측 모델 개발", 한국전자파학회지 전자파기술, 제10권, 제1호, pp. 2–15, 1999년 3월.
[6] 손해원, 도심 마이크로셀과 실내 전파 환경 예측을 위한 결정적인 레이 튜브 방법 (A Deterministic Ray Tube Method for Wave Propagation Predictions in Urban Microcellular and Indoor Environments), KAIST 박사 학위 논문, 2000. (방문일 2022-08-20)
[7] 이행선, "실외 전파 특성 계산을 위한 고속 3차원 광선 추적법", 한국전자파학회논문지, 제18권, 제11호, pp. 1231–1236, 2007년 11월.
[8] H. Kim and H.-S. Lee, "Accelerated three dimensional ray tracing techniques using ray frustums for wireless propagation models," Prog. Electromagn. Res., vol. 96, pp. 21–36, 2009.
[9] 권세웅, 문현욱, 오재림, 임재우, 배석희, 김영규, 박정수, 윤영중, "가속 방법을 이용하는 전파 광선 추적법에 관한 연구", 한국전자파학회논문지, 제20권, 제5호, pp. 471–479, 2009년 5월.
[10] O. Landron, M. J. Feuerstein and T. S. Rappaport, "A comparison of theoretical and empirical reflection coefficients for typical exterior wall surfaces in a mobile radio environment," IEEE Trans. Antennas Propag., vol. 44, no. 3, pp. 341–351, Mar. 1996.
[11] H. M. El-Sallabi and P. Vainikainen, "Improvements to diffraction coefficient for non-perfectly conducting wedges," IEEE Trans. Antennas Propag., vol. 53, no. 9, pp. 3105–3109, Sep. 2005.
[12] 윤대영, 이해승, 이소영, 박용배, "LBVH 순회를 이용한 거대 도체 구조의 레이다 단면적 분석", 한국전자파학회논문지, 제33권, 제9호, pp. 739–742, 2022년 9월.
[13] G. Carluccio and M. Albani, "An efficient ray tracing algorithm for multiple straight wedge diffraction," IEEE Trans. Antennas Propag., vol. 56, no. 11, pp. 3534–3542, Nov. 2008.
[14] V.-D. Nguyen, A. Mansour, A. Coatanhay, and T. Marsault, "A generalized recursive Vogler algorithm for multiple bridged knife-edge diffraction," IEEE Trans. Antennas Propag., vol. 70, no. 11, pp. 10828–10842, Nov. 2022.

[다음 읽을거리]

2022년 8월 8일 월요일

반복적 물리 광학(反復物理光學, IPO: Iterative Physical Optics)

[경고] 아래 글을 읽지 않고 "반복적 물리 광학"을 보면 바보로 느껴질 수 있습니다.

[확인] 본 페이지는 exp(-iωt) 시간 약속을 사용하고 있습니다.


물리 광학(物理光學, physical optics, PO) 혹은 PO가 다루지 못하는 다중 반사(multiple reflection)를 어느 정도 정확하게 계산하기 위해 제안된 기법은 반복적 물리 광학(反復物理光學, iterative physical optics, IPO)[1] 혹은 IPO이다. 물리 광학은 전자파가 조사된 PEC 표면만 고려해서 산란을 계산하므로, 다중 반사가 일어나는 상황에서는 오차가 심해진다. 반면에 반복적 물리 광학은 물리 광학으로 계산을 시작하지만, 각 표면이 가진 전류 밀도로 산란장을 만들어서 다시 기존의 전류 밀도를 갱신하는 방식으로 전류 밀도의 오차를 줄여나간다. 이 과정은 사용자가 원하는 만큼 반복할 수 있어서 반복적 물리 광학으로 부른다.
반복적 물리 광학의 공식화는 스트래튼–추 공식(Stratton–Chu Formula)으로 만든 MFIE(자기장 적분 방정식, Magnetic Field Integral Equation)부터 출발한다.

                        (1)

여기서 $g(\bar r, \bar r')$는 3차원 자유 공간 그린 함수(3D free-space Green's function) $g(\bar r, \bar r'; k)$와 동일하다. 식 (1)에서 $\bar r$ = $\bar r'$인 경우가 있으므로, 식 (1)에 나온 표면 적분은 $\bar r$ = $\bar r'$에서 [그림 1]과 같은 특이점을 항상 가지고 있다.

[그림 1] 미소 면적소 $s_0'$에 대한 좌표계

따라서 식 (1)을 제대로 계산하려면, 특이점이 있는 미소 면적소 $s_0'$의 포함 여부에 따라 적분을 분리해서 연산해야 한다.

                        (2)

면적 $s' - s_0'$에는 $\bar r$ = $\bar r'$인 특이점이 없어서 기존 적분 방식을 쓰면 된다. 수학적으로 고민되는 부분은 [그림 1]에서 $R \to 0$으로 가는 미소 면적소 $s_0'$에 대한 적분이다. 이 적분을 위해 $\bar J_s(\bar r)$을 접선 자기장 $\bar H_t (\bar r)$로 분해한 후 벡터 삼중적(vector triple product)을 적용해서 정리한다.

                  (3a)

                        (3b)

여기서 $\bar r - \bar r'$ = $R \hat R$, $R \to 0$, $\hat n'$과 $\bar \nabla' g(\bar r, \bar r')$은 항상 수직이다. 식 (3)을 식 (2)에 대입해서 $\bar J_s (\bar r)$에 대해 정리하면, PEC에 사용되는 반복적 물리 광학을 위한 표현식이 얻어진다.

                        (4)

그린 함수 $g (\bar r, \bar r')$의 구배(gradient)를 구하기 위해, 구의 중심을 $\bar r'$으로 선택해서 편미분한다.

                        (5a)

                        (5b)

전류 밀도 $\bar J_s (\bar r)$이 반복될 수 있도록 식 (4)를 바꾼다[1].

                        (6)

여기서 $m$ = $0, 1, 2, \cdots$, $\bar J_0 (\bar r)$ = $2 \hat n \times \bar H_i(\bar r)$이다. 물리 광학 혹은 PO처럼 구조물을 삼각화(三角化, triangulation)한 후에, 식 (6)에 따라 반복적 물리 광학 혹은 IPO를 적용하는 순서는 다음과 같다.
  • 모든 삼각형의 전류 밀도를 $0$으로 만든다.
  • 전자파에 조사된 삼각형의 전류 밀도는 $\bar J_0 (\bar r)$ = $2 \hat n \times \bar H_i(\bar r)$로 초기화한다.
  • 크기가 $0$이 아닌 전류 밀도 $\bar J_m (\bar r)$에 식 (6)의 우변을 적분해서 $\bar J_{m+1} (\bar r)$을 얻는다. 다만 $\bar r'$에서 $\bar r$로 가는 광선은 가려짐 없이 $\bar r$이 가르키는 삼각형을 조사해야 한다.
  • 원하는 정밀도가 나올 때까지 $\bar J_{m+1} (\bar r)$ 계산을 반복한다. 
IPO는 역행렬을 계산하는 과정 없이 반복으로만 답을 구하기 때문에, 복잡한 구조물의 산란 문제를 탐구할 때에 적절한 기법이다. 하지만 구조물의 크기가 매우 커져서 더해지는 항이 많아지면, 행렬의 반복법(iterative method of matrix)에서 경험한 현상처럼 $\bar J_{m+1} (\bar r)$이 발산하는 문제가 필연적으로 생긴다.
임피던스 경계 조건(impedance boundary condition)을 써서 IPO의 적용 범위를 더 다양한 매질로 확대할 수 있다[3], [4]. IPO 공식화는 임피던스 경계 조건에 대한 MFIE(자기장 적분 방정식, Magnetic Field Integral Equation)에서 시작한다.

                  (7)

여기서 임피던스 경계면의 등가 자류 밀도는 $\bar M_s(\bar r')$ = $Z_s \bar J_s(\bar r') \times \hat n'$이다. 식 (7)에 식 (3b), (5b)를 적용해서 적분 방정식을 풀 수 있는 형태로 바꾼다. 다만 미분 연산자 $\bar M_s(\bar r') \cdot \bar \nabla'$은 공식화가 쉬운 데카르트 좌표계(Cartesian coordinate system)에서 풀어서 정리한다.

                  (8a)

                  (8b)

                  (9)

여기서 $g$ = $g(\bar r, \bar r')$, $\bar M_s$ = $\bar M_s(\bar r')$ = $M_x \hat x + M_y \hat y + M_z \hat z$이다. 식 (8b)에 나온 그린 함수의 2번 미분 $\partial^2 g(\bar r, \bar r')/ \partial R^2$은 다음과 같이 계산된다.

                  (10)

점 $\bar r$ = $\bar r'$ 근방 혹은 $s_0'$ 표면에서, 식 (9)의 둘째줄과 셋째줄에 있는 $Z_s$가 곱해진 표면 적분은 식 (3)과 다르게 $0$으로 간다. 왜냐하면 $R \to 0$일 때, $g(\bar r, \bar r')$의 크기는 $\partial g(\bar r, \bar r') \mathbin{/} \partial R$보다 천천히 발산하고, $\bar M_s(\bar r')$ 혹은 $\bar J_s(\bar r') \times \hat n'$의 발산은 전류 및 자류 밀도가 표면에만 있어서 $0$이 되기 때문이다. 식 (6)과 비슷하게 식 (9)를 IPO를 위한 재귀 방정식으로 변환한다[4].

                  (11)

여기서 $m$ = $0, 1, 2, \cdots$이다.

[참고문헌]
[1] F. Obelleiro-Basteiro, J. L. Rodriguez and R. J. Burkholder, "An iterative physical optics approach for analyzing the electromagnetic scattering by large open-ended cavities," IEEE Trans. Antennas Propag., vol. 43, no. 4, pp. 356–361, Apr. 1995.
[2] 이용희, 진희철, 김경태, "GPU를 이용한 반복적 물리 광학법의 가속화에 대한 연구", 한국전자파학회논문지, 제26권, 제11호, pp. 1012–1019, 2015년 11월.
[3] A. W. Glisson, "Electromagnetic scattering by arbitrarily shaped surfaces with impedance boundary conditions," Radio Sci., vol. 27, no. 06, pp. 935–943, Nov.–Dec. 1992.
[4] J.-W. Rim and I.-S. Koh, "Convergence and accuracy of near-field-corrected iterative physical optics for scattering by imperfectly conducting and dielectric objects," IET Microw. Antennas Propag., vol. 14, no. 10, pp. 999–1005, Aug. 2020.
[5] M. A. Shah, Ç. Tokgöz, and B. A. Salau, "Radar cross section prediction using iterative physical optics with physical theory of diffraction," IEEE Trans. Antennas Propag., vol. 70, no. 6, pp. 4683–4690, Jun. 2022.

2022년 8월 3일 수요일

물리 광학(物理光學, PO: Physical Optics)

[경고] 아래 글을 읽지 않고 "물리 광학"을 보면 바보로 느껴질 수 있습니다.

[확인] 본 페이지는 exp(-iωt) 시간 약속을 사용하고 있습니다.


물리 광학(物理光學, physical optics, PO) 혹은 PO는 고전적 기하 광학(geometrical optics)빛이 편광(偏光, polarization)을 가진 파동이란 개념을 더해서 전자파의 산란(散亂, scattering)을 더 정확하게 계산하는 방법이다. 계산의 정밀도 관점으로 보면, 기하 광학 $<$ 물리 광학 $<$ 맥스웰 방정식(Maxwell's equations) 순이다. 완전 정확한 맥스웰 방정식의 해법이 있는 현시점에 근사가 심한 기하 광학이나 물리 광학이 필요할까? 현업에서는 기하 광학과 물리 광학의 존재 가치가 크다. 작은 문제에 대해 우리가 원하는 수준까지 정확하게 맥스웰 방정식을 풀 수 있지만, 모든 문제를 맥스웰 방정식에 대입해서 해결할 수는 없다. 대개 파장에 비해 매우 큰 산란체는 맥스웰 방정식으로 해결할 수 없는 경우가 대부분이다. 왜냐하면 구조물을 이산화할 때 생기는 미지수는 3차원 부피 문제에서 $n^3$ 비율로 커지기 때문이다. 여기서 $n$은 한 축당[$x, y, z$중 하나] 이산화에 사용하는 기저(basis)의 개수이다. 경험 법칙(rule of thumb) 관점에서 산란체의 축당 이산화 개수는 대략 $n$ $\approx$ $D / (\lambda/10)$이다. 여기서 $D$는 산란체의 크기, $\lambda$는 매질속의 파장이다. 또한 미지수로 연립 방정식(simultaneous equations)을 구성해서 풀 때 필요한 행렬의 차원(dimension)은 보통 $n^3 \times n^3$이다. 그러면 행렬을 저장하기 위한 공간만 $n^6$이 필요하다. 축당 기저 개수 $n$이 작을 때는 느끼기 어렵지만, 구조물이 파장보다 매우 커지면 저장 공간과 계산 시간이 기하 급수적으로 커진다. 예를 들어, 산란체의 길이 $D$가 $10 \lambda$라면, $n$ = $100$이 되어서 저장 공간은 $10^{12}$개가 필요하다. 계산에 배정밀도 부동 소수점 형식(double-precision floating-point format)을 선택한다면, 공간 하나당 8바이트(byte)가 필요해서 최종 저장 공간은 행렬 하나에 8 TB가 된다. 3차원 부피가 아닌 2차원 표면만 고려하면, 저장 공간은 $n^2 \cdot n^2$ = $n^4$이라서 최종 8 MB만 필요하다. 하지만 $D$ = $100 \lambda$이면, 저장 공간은 다시 8 TB로 확 늘어난다. 따라서 구조물이 매우 큰 경우에는 맥스웰 방정식을 직접 풀기가 불가능해서 기하 광학이나 물리 광학의 적용이 필수적이다.

[그림 1] 구조물에 의해 산란되는 전자기장

고전적(classical) 물리 광학은 파동 개념을 쓰기 때문에 기하 광학보다는 물리에 더 가까운 광 이론이란 개념으로 쓰였다. 하지만 현재는 맥스웰 방정식이란 정답이 있으므로 기존보다  적용 범위를 더 좁혀서, 양자 역학(quantum mechanics)에 쓰는 보른 근사(Born approximation)[1] 혹은 전자파를 위한 키르히호프 근사(Kirchhoff approximation)를 사용하여 산란파를 대충 계산하는 방식을 물리 광학 혹은 PO라 부른다. 보른 근사는 양자 상태(quantum state)를 나타내는 파동 함수(wave function)의 산란파 $U_s$를 정확히 계산하지 않고 입사파 $U_i$와 산란파의 합인 전체파 $U_\text{tot}$를 입사파로 대충 근사화하는 방식이다.[$U_\text{tot}$ = $U_i + U_s$ $\approx$ $U_i$] 보른 근사와 비슷하게, 전자파 산란이 일어나는 경우에 키르히호프 근사는 산란파의 접선 및 법선 성분을 입사파로 바꾸어서 계산한다. 키르히호프 근사를 위해 선택하는 산란파와 입사파의 경계 조건은 키르히호프 경계 조건(Kirchhoff boundary condition)이라 부른다. 대부분의 물리 광학 문제는 완전 자기 도체(Perfect Electric Conductor, PEC) 상에 유기되는 표면 전류 밀도 $\bar J_s(\bar r)$을 입사 자기장 $\bar H_i(\bar r)$로 근사한다. 즉, PEC에서 키르히호프 경계 조건은 $\hat n \times \bar H_s(\bar r)$ $\approx$ $\hat n \times \bar H_i(\bar r)$이다. 보른 근사 혹은 키르히호프 근사처럼 물리 광학은 연립 방정식을 풀지 않기 때문에 파동의 산란 현상을 매우 빠르게 예측하지만 약간 오차가 있다[2], [3]. 물리 광학의 밑바탕을 이해하기 위해, [그림 1]과 같은 상황에서 PEC에 의해 생기는 산란 자기장 $\bar H_s(\bar r)$을 정확하게 기술한다.

                   (1)

                   (2)

                   (3)

여기서 $s'$는 모든 PEC 표면적, $\hat n$은 PEC 표면에 수직인 단위 벡터, 자기 벡터 포텐셜 $\bar A(\bar r)$에 대한 3차원 자유 공간 그린 함수(3D free-space Green's function)는 $G_A(\bar r, \bar r'; k)$ = $e^{i k R} \mathbin{/} (4 \pi R)$, $k$ = $\omega \sqrt{\mu \epsilon}$, $R$ = $|\bar r - \bar r'|$, $\bar r - \bar r'$ = $R \hat R$이다. 원래 $\bar J_s(\bar r)$는 입사와 산란 자기장에 의해 생기지만, 키르히호프 근사에 따라 산란 자기장은 거의 입사 자기장과 같은 $\bar H_s(\bar r)$ $\approx$ $\bar H_i(\bar r)$이라고 가정한다. 이는 기하 광학의 광선이 산란 없이 직선으로 입사하여 PEC에 만드는 전류 밀도와 동일해서 기하 광학 전류 밀도(GO electric current density) $\bar J_\text{GO}(\bar r)$이라고 부른다. 

                   (4)

여기서 $2\bar H_i(\bar r)$ $\approx$ $\bar H_i(\bar r) + \bar H_s(\bar r)$이다. 식 (4)를 식 (3)에 대입해서 다시 정리한다.

             (5)

여기서 $\partial / \partial n'$ = $\hat n' \cdot \bar \nabla'$이다. 빛이 조사되어 물리 광학에 기여하는 PEC 표면 $a_i'$은 $\hat s \cdot \hat n' < 0$ 조건으로 선택한다.[단위 계단 함수로 쓰면 $u(-\hat s \cdot \hat n')$] 여기서 $\hat s$는 전자파의 전력 전달 방향을 나타내는 포인팅 단위 벡터(Poynting unit vector)이다. 근사적으로 유도한 식 (5)조차 복잡하므로, $|\bar r| \gg |\bar r'|$인 원역장 조건(far-field condition)을 써서 더욱 간략화한다.

                   (6)

여기서 $\bar r$ = $(r, \theta, \phi)$, $\bar \nabla$ = $i \bar k$ = $ik \hat r$, $\eta$는 매질의 고유 임피던스(intrinsic impedance)이다. 그린 함수 $G_A(\bar r, \bar r'; k)$에도 원역장 근사(far-field approximation)를 적용한다.

                  (7)

여기서 $\psi (\bar r')$은 광선이 PEC 위치 $\bar r'$에 도달할 때 생기는 원점 대비 위상 차이이다.

[그림 2] 삼각화를 이용한 그물망 생성(출처: wikipedia.org)

[그림 3] 제$m$번 삼각형의 구조

[그림 2]처럼 산란체를 삼각화하여 표면 적분을 유한 합으로 근사함으로써 식 (7)의 마지막식을 더 간단하게 나타낸다.

                   (8)

여기서 $M$은 물체를 구성하는 그물눈(mesh)의 개수, $u(\cdot)$는 단위 계단 함수(unit step function), $\bar c_m$과 $\Delta a_m$ 각각은 [그림 3]에 있는 제$m$번 그물눈을 구성하는 삼각형의 무게 중심(barycenter)과 면적, $\hat n$은 이 삼각형 면적의 단위 벡터, 삼각형의 면적 벡터는 $\Delta \bar a_m$ = $\Delta a_m \hat n$이다. [그림 2]에 보인 삼각화(三角化, triangulation)는 구조체를 형상화하기 위한 그물망 생성(mesh generation)의 대표적인 방법이다. 물리 광학 계산이 정밀도를 유지하려면, 경험적으로 [그림 3]에 그린 삼각형의 한 변 길이가 $\lambda/10$보다 작아야 한다. [그림 3]에서 제$m$번 삼각형의 무게 중심 $\bar c_m$과 면적 $\Delta a_m$은 좌표계 기반 벡터(vector)로 쉽게 구할 수 있다.

                   (9)

또한 함수 $u(\hat r \cdot \hat n)$은 표면적의 바깥 방향인 $\hat n$ 성분을 가지고 $\hat r$방향으로 전파하는 자기장만 물리 광학의 산란에 기여한다는 뜻이다. 다만 전자파에 조사된 표면적을 선택하는 $u(-\hat s \cdot \hat n')$ 조건과 헷갈리면 안된다. 단위 벡터 $\hat s$는 입사파가 들어오는 방향, $\hat r$은 산란파가 진행하는 방향이다. 즉, 표면적 $\Delta a_m$에 들어오고 나가는 전자파의 방향이 다르기 때문에, $u(-\hat s \cdot \hat n')$과 $u(\hat r \cdot \hat n)$처럼 부호를 바꾸어서 선택한다. 만약 $\bar r$이 $\bar r'$보다 아주 크지 않고 $\bar r'$의 영향을 고려해야 한다면, 식 (8)의 $\bar H_s(\bar r)$과 $\Delta \bar P_m(\bar r)$은 다음처럼 바뀐다.

                   (10)

여기서 $\hat R$은 $\bar r - \bar r'$의 단위 벡터이다.

[그림 4] 산란체를 표면 전류 및 자류 밀도로 등가화

PEC라고 생각한 산란체가 임의의 매질로 일반화된다면, $\mu$와 $\epsilon$이 마음대로 변해도 반사를 계산하는 프레넬 방정식(Fresnel equation)이 필요하다. 키르히호프 근사와 [그림 4]에 소개한 표면 등가의 원리(surface equivalence principle)에 따라 기하 광학 전류 밀도 $\bar J_\text{GO}(\bar r)$과 자류 밀도 $\bar M_\text{GO}(\bar r)$를 근사화한다. 입사파가 표면에 거의 수직하게 들어오면, $\bar J_\text{GO}(\bar r)$과 $\bar M_\text{GO}(\bar r)$은 매우 간단하게 구해진다.

                   (11)

여기서 표면에 대한 입사각은 $\theta_i \approx 0$, $r_p$와 $r_s$는 각각 P편광[혹은 평행 편파, TM파]S편광[혹은 직각 편파, TE파]반사 계수(reflection coefficient)이다. PEC 경우에 $\epsilon_\text{in} \to \infty$로 가서 $r_p$ = $1$, $r_s$ = $-1$이므로, 식 (11)은 식 (4)에 수렴한다. 따라서 전기 원천 $\bar J_\text{GO}(\bar r)$과 자기 원천 $\bar M_\text{GO}(\bar r)$이 만드는 전자기장은 식 (8)과 비슷하게 표현된다.

                   (12)

                   (13)

                   (14)

여기서 $\epsilon_\text{in}$은 산란체 내부의 유전율, $G_F(\bar r, \bar r'; k)$ = $G_A(\bar r, \bar r'; k)$, $(\mathsf{e})$와 $(\mathsf{m})$은 각각 전기와 자기 원천을 의미한다. 맥스웰 방정식의 쌍대성(duality)을 적용해서 식 (12)를 식 (13)으로 손쉽게 바꿀 수도 있다. 최종적으로 임의 매질에 의해 산란되는 전자기장은 다음과 같다.

                   (15)

물리 광학이 좋은 근사이기는 하지만, 기본적으로 기하 광학과 같은 근축 근사(paraxial approximation)이다. 그래서 반사면에 수직이어서 반사가 많이 일어나는 방향에서는 물리 광학이 잘 맞지만, 반사면에 평행하거나 반대편으로 생기는 산란이나 회절(回折, diffraction)은 물리 광학으로 예측할 수 없다. 수직 입사(normal incidence)에서 많이 벗어나는 경우는 광선 고정 좌표계(ray-fixed coordinate system)를 써서 식 (11)을 더 정확하게 공식화한다.

                   (16)

여기서 $\hat e_\parallel^i$와 $\hat e_\perp$는 입사 평면(plane of incidence)에 평행 및 수직이고 입사 광선(incident ray)에는 수직인 단위 벡터, $\hat e_\parallel^r$은 반사 평면(plane of reflection)에 평행하면서 반사 광선(reflection ray)에는 수직이다.

[참고문헌]
[1] M. Born, "Quantenmechanik der Stoßvorgänge (Quantum mechanics of collision processes)," Zeitschrift für Physik (Magazine for Physics), vol. 38, pp. 803–827, Nov. 1926.
[2] K. M. Siegel, H. A. Alperin, R. R. Bonkowski, J. W. Crispin, A. L. Maffett, C. E. Schensted, and I. V. Schensted, "Bistatic radar cross sections of surfaces of revolution, J. Appl. Phys., vol. 26, no. 3, 297–305, Mar. 1955.
[3] R. G. Kouyoumjian, "Asymptotic high-frequency methods," Proc. IEEE, vol. 53, no. 8, pp. 864–876, Aug. 1965.
[4] I. Gupta and W. Burnside, "A physical optics correction for backscattering from curved surfaces," IEEE Trans. Antennas Propag., vol. 35, no. 5, pp. 553–561, May 1987.
[5] R. Cicchetti and A. Faraone, "Analysis of open-ended circular waveguides using physical optics and incomplete Hankel functions formulation," IEEE Trans. Antennas Propag., vol. 55, no. 6, pp. 1887–1892, Jun. 2007.
[6] X. F. Li, Y. J. Xie, P. Wang and T. M. Yang, "High-frequency method for scattering from electrically large conductive targets in half-space," IEEE Antennas Wirel. Propag. Lett., vol. 6, pp. 259–262, 2007.
[7] 석성하, 물리 광학 방법을 이용한 레이더 단면적 예측에 관한 연구 (Estimation of Radar Cross Section Using Physical Optics), POSTECH 석사 학위 논문, 1997.
[8] A. W. Glisson, "Electromagnetic scattering by arbitrarily shaped surfaces with impedance boundary conditions," Radio Sci., vol. 27, no. 06, pp. 935–943, Nov.–Dec. 1992.
[9] 신호근, 박용배, "물리광학법과 물리광학 회절이론을 이용한 레이다 단면적 해석 기법", 한국전자파학회지 전자파기술, 제28권, 제6호, pp. 34–39, 2017년 11월.
[10] F. L. Song, Y. M. Wu, Y. Pan, J. Hu, and Y. -Q. Jin, "Fast physical optics method based on adaptive mesh technique," IEEE Trans. Antennas Propag., vol. 71, no. 11, pp. 9113–9118, Nov. 2023.

[다음 읽을거리]