2022년 7월 5일 화요일

장거리 전자파 전파 모형화(Long-range Electromagnetic Propagation Modeling)

[경고] 아래 글을 읽지 않고 "장거리 전자파 전파 모형화"를 보면 바보로 느껴질 수 있습니다.

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


맥스웰 방정식(Maxwell's equations)은 고전적 의미에서 완벽하기 때문에 원론적으로 풀지 못하는 문제는 없다. 하지만 현실에서는 정밀도와 계산 시간 한계로 인해 맥스웰 방정식만으로 문제를 해결할 수는 없다. 특히 원천에서 수신부가 많이 멀어지는 경우에는 편미분 방정식(partial differential equation, PDE)인 맥스웰 방정식을 있는 그대로 적용하지 않고, 물리학에 바탕을 둔 광선 추적(ray tracing)을 사용해서 반복법 형태로 답을 구한다.

[그림 1] 지구 대기권의 특성(출처: wikipedia.org)

광선 추적도 통하지 않는 대류권(對流圈, troposphere) 영역과 같은 진짜 장거리 문제를 풀려면, 맥스웰 방정식을 적절히 근사해서 정확도와 계산 시간을 적절히 타협한 현실적인 방정식 모양으로 바꾸어야 한다[1]. 장거리에서 사용 가능한 전자파 전파를 모형화하기 위해 맥스웰 방정식부터 출발한다. 대기를 다루고 있어서 유전율은 $\epsilon$ = $\epsilon(\bar r)$과 같이 공간적으로 변하고 투자율은 $\mu$ = $\mu_0$로 일정하다.

                  (1a)

                  (1b)

                  (2)

여기서 $k$ = $\omega \sqrt{\mu \epsilon}$이다. 식 (1), (2)에 유도한 파동 방정식을 간략화하기 위해, 방위각 $\phi$방향으로는 전자기장이 일정하다고 가정하고[$\partial / \partial \phi$ = $0$], 3차원 공간 전체를 효율적으로 다루는 구 좌표계(spherical coordinate system)도 도입한다. 그러면 $E_\phi$는 $H_r, H_\theta$, $H_\phi$는 $E_r, E_\theta$를 각각 생성하므로, 전자기장은 $E_\phi, H_\phi$를 기준으로 완벽히 분해된다. 먼저 어려워보이는 식 (1b)에 유일한 자기장 성분인 $H_\phi$를 넣고 정리한다[1].

                  (3)

여기서 $\partial_r$ = $\partial / \partial r$, $\partial_\theta$ = $\partial / \partial \theta$, $\epsilon$ = $\epsilon(r, \theta)$이다. 지면을 $xy$평면으로 가정한 경우, 식 (3)은 VED(수직 전기 쌍극자, vertical electric dipole)가 만드는 수직 편파(vertical polarization, VP)를 뜻한다. 식 (2)는 식 (1)의 쌍대(duality)라서, 식 (3)에 따라 전기장 성분 $E_\phi$의 방정식을 얻는다.

                  (4)

전기장 $E_\phi$는 $xy$평면에 평행인 전류 고리가 만들기 때문에, 식 (4)는 HED(수평 전기 쌍극자, horizontal electric dipole)에 의한 수평 편파(horizontal polarization, HP)를 나타낸다. 쌍대를 써서 HED 대신 VMD(수직 자기 쌍극자, vertical magnetic dipole)라 표현하기도 한다. 최종적으로 식 (3)과 (4)가 장거리 전자파 전파를 위한 편미분 방정식이 된다.
공간에서 유전율이 $\epsilon$ = $\epsilon(r, \theta)$로 변함으로 인해 식 (3)의 해를 완벽히 구하기는 불가능하다. 그래서 $r \gg 1$인 조건에서 2차원과 3차원에 대한 급속 하강 방법(method of steepest descent)의 결과를 이용한다.

                      (5)

                      (6)

식 (5)와 (6)을 참고해서 $H_\phi$를 대략적으로 공식화한다.

                  (7)

왜냐하면 식 (5)에 따라 $H_\phi$는 2차원 기준으로 $1/\sqrt{\rho}$에 비례하고 식 (6)에 나온 3차원 기준으로 $k$에도 비례하기 때문이다. 여기서 $\rho$ = $r \sin \theta$, $U(r, \theta)$는 2차원만 고려한 파동이다. 조금 더 설명하면, 방위각 대칭성을 사용해서 식 (6)의 적분 변수를 다음처럼 살짝 바꾼다.

                      (8)

그러면 $g(\zeta, \xi)$를 $\zeta$에 대해 적분함으로써 $1/\sqrt{\rho}$가 나오고 $U(r, \theta)$는 $\xi$에 대한 푸리에 변환으로 나온 결과에 해당한다. 다만 $k \sin \theta \sin \phi$ = $\bar k \cdot \hat y$는 식 (8)의 최종 항이 $y$방향에 평행임을 표현하고 있어서 $1/\sqrt{\sin \theta}$를 상쇄시키지 못한다.
편미분 방정식의 관점에서 식 (8)의 의미를 파악하기 위해 $H_\phi$을 일반식으로 가정해 식 (3)에 대입한다.

                      (9)

여기서 $\alpha, \beta, \gamma$는 적당한 상수이다. 그러면 식 (3)의 첫째 항은 다음 관계를 만족한다.

                      (10a)

                      (10b)

만약 식 (7)처럼 $\alpha$ = $1/2$라면, 식 (10b)는 통상적인 헬름홀츠 방정식(Helmholtz equation) 형태로 바뀐다.

                      (11a)

                      (11b)

식 (3)의 둘째 항은 다소 복잡해서 식 (10b)처럼 $\partial_\theta \epsilon \partial_\theta U$를 가진 항부터 본다. 우선 식 (3)의 둘째 항에 나온 $\partial_\theta$부터 계산한다.

                      (12)

식 (12)에 다시 $\partial_\theta$를 적용해서 $\partial_\theta \epsilon \partial_\theta U$ 항만 모은다.

                      (13)

더 나가서 식 (3)의 둘째 항을 $\theta$에 대한 헬름홀츠 방정식으로 만들기 위해, $\partial_\theta U$ 항이 0이 되는 조건도 찾는다.

                      (14)

식 (13)과 (14)에 찾은 $\alpha$ = $1/2$, $\beta$ = $-1/2$ 조건을 식 (12)에 넣고 $\theta$에 대해 편미분한다.

                      (15a)

                      (15b)

식 (11b)와 (15b)를 식 (3)에 넣어서 $U(r, \theta)$에 대한 새로운 헬름홀츠 방정식을 유도한다[1].

                      (16a)

                      (16b)

상수 $\gamma$의 선택에는 특별한 제한이 없지만 식 (5)에 따라 $\gamma$ = $-1/2$로 택한다. 최종적으로 수직 편파 혹은 VED에 대한 장거리 전자파 전파 모형화(long-range electromagnetic propagation modeling)가 완성된다.

                  (17)

여기서 $k_0$ = $\omega \sqrt{\mu_0 \epsilon_0}$, 굴절률(refractive index)은 $n(r, \theta)$ = $\sqrt{\epsilon(r, \theta) / \epsilon_0}$이다. 또한 식 (7)을 참고해서 전기장 $E_\phi$의 공식은 다음처럼 선택한다.

                  (18)

식 (18)에 따라 수평 편파 혹은 HED를 위한 장거리 전자파 전파 모형화도 얻는다.

                  (19)

[그림 2] 방위각 $\phi$ = $0^\circ$인 구 좌표계 $(r, \theta)$를 데카르트 좌표계 $(x, z)$로 변환

구 좌표계가 전자파 전파 모형화에 유리하지만 편미분 계산이 어려워지고 지표면 평평하다는 우리의 직관과 맞지 않아서, [그림 2]처럼 방위각을 $\phi$ = $0^\circ$[방위각 기준점을 편하게 0˚로 설정]으로 둔 구 좌표계 $(r, \theta)$를 평평한 지구로 근사화한 데카르트 좌표계 $(x, z)$로 바꾼다[1].

                      (20)

여기서 $R_e$는 지구의 평균 반지름이며 대략 6,371 km이다. 식 (20)을 써도 되지만, $r$이 고정되고 $\theta$를 변화시킬 때는 [그림 2]처럼 수평 방향 $x$를 기준으로 높이 $z$는 지구의 곡률만큼 다소 크게 변한다. 이를 완화시키기 위해 뫼비우스 변환(Möbius transformation)을 도입한다. 스미쓰 도표(Smith chart)에도 쓰이는 뫼비우스 변환은 복소 평면에서 원의 내부와 무한 반평면을 등각 사상(conformal mapping)으로 연결한다. 그래서 식 (20)을 변형해서 뫼비우스 변환 모양으로 바꾸면, $x$에 대해 $z$가 휘어지는 정도를 감소시킬 수 있다.

                  (21)

여기서 $\xi$ = $x + iz$, $\gamma$ = $r \sin \theta + i r \cos \theta$이다. 식 (21)은 전형적 뫼비우스 변환인 식 (21)은 복소 변수 $\gamma$를 기준으로 판별식 $ad-bc$ = $-2R_e$ $\ne$ $0$이므로 절대 약분되어 상수가 될 수 없다. 식 (21)을 써서 식 (17)을 좌표 변환하기 위해 $(r, \theta)$를 위한 2차원 라플라시안 $\nabla_{r \theta}^2$을 사용한다.

                      (22a)

                      (22b)

                      (22c)

식 (22b)를 편미분해서 더 사용하기 편한 형태로 바꾼다.

                      (23)

식 (21)에 따라 식 (22c)를 $(x, z)$ 좌표계로 다시 쓴다.

                      (24a)

                      (24b)

전자파가 거의 $x$축을 따라서 전파된다고 생각해서 $U(x, z)$ = $e^{i k_0 x} u(x, z)$로 치환한다[1].

                  (25)

여기서 $n_d$는 굴절률 $n$의 등가적 변화량이다. 편미분 방정식의 분류 기준에서 식 (25)는 포물형 편미분 방정식(parabolic partial differential equation)이라서 간단히 포물형 방정식(parabolic equation, PE)으로 부르기도 한다. 왜냐하면 원뿔 곡선의 판별식(discriminant of conic section) $D$ = $ac - b^2$ = $1\cdot 1 - 1^2$ = $0$이 되기 때문이다. 또한 식 (25)는 근축 근사(paraxial approximation)의 결과와도 연결되어 근축 파동 방정식(paraxial wave equation)의 범주에 들어간다. 비슷한 방식으로 식 (19)도 데카르트 좌표계에서 공식화할 수 있다.

                  (26a)

                      (26b)

여기서 $V(x, z)$ = $e^{i k_0 x} v(x, z)$이다. 구와 데카르트 좌표는 식 (21)을 정리하고 실수부와 허수부를 등식화해서 구한다.

                      (27)

여기서 $R^2$ = $r^2 + R_e^2 + 2 R_e r \cos \theta$이다. 식 (27)은 다소 복잡해서 $\theta \approx 0$이라는 조건으로 근사화한다.

                      (28)

여기서 $R \approx r + R_e$, $dx \approx r d \theta$, $dz \approx dr$이다. 굴절률 $n$의 공간 변화는 거의 0이라 두고 식 (28)을 식 (25)에 대입한다.

                      (29a)

                      (29b)

식 (29a)에 나온 $m_d$도 $n_d$와 비슷하기 때문에[$m_d \approx n_d$], 전자파의 장거리 전파 특성에는 편파 영향이 거의 없다. 식 (29b)는 지구를 평평하다고 가정해서 나온 항 $2z/R_e$를 포함하고 있어서, 이 항을 강조한 등가 굴절율 $n_\text{eff}$를 다음과 같이 정의할 수 있다.

                      (30a)

                      (30b)

                  (30c)

여기서 $n$ = $1 + 10^{-6}N$, $N$과 $M$은 각각 굴절도(refractivity)변형 굴절도(modified refractivity), $h$는 지표면에서 계산 지점까지의 높이, $10^6$은 대기 굴절도를 보기 좋게 만드는 상수이다. 지구는 실제 둥글지만 평평하다고 강제로 가정하면, 전방으로 쏜 전자파가 마치 하늘로 올라가는 것처럼 보인다. 이러한 왜곡을 굴절율로 보정하는 항이 변형 굴절도 $M$에 있는 $h/R_e$이다. 등각 사상 관점에서 $h/R_e$은 구 좌표계에서 데카르트 좌표계로 갈 때에 생기는 길이 증가분이다. 하지만 이 개념은 어려우니까 전자파가 평평한 지면을 따라 진행할 때 하늘로 굴절되는 양으로서 $h/R_e$를 굴절도에 추가한다. 굴절도 $N$을 무시하면 실제로 전자파는 직선 방향으로 잘 진행하며, 둥근 지구로 인해 전자파가 아닌 땅이 휘어진다. 대기중의 전파 특성 분석에는 높이에 대한 굴절도 기울기도 많이 쓴다. 보통 $dN/dh$, $dM/dh$를 각각 N구배(N-gradient)M구배(M-gradient)라고 이름 붙인다. 높이 올라갈수록 대기가 약해져서 N구배는 음수이며 표준값은 대략 km당 -39[= -39 N/km] 정도가 된다. 물론 지역마다 대기 특성이 크게 다르기 때문에 N구배는 지역, 계절, 높이 등에 따라 심하게 변한다. M구배는 식 (30c)에 따라 $dN/dh + 10^6/R_e$로 정확히 정의된다. 여기서 $10^6/R_e \approx 157$ M/km이다. N구배의 표준값으로 환산한 M구배는 대충 118 M/km가 된다. N구배는 음수라서 전자파는 땅으로 구부러지려 하지만, $10^6/R_e$이 매우 크기 때문에 최종적인 M구배는 양수가 되며 전자파는 하늘로 약간 굽어지게 된다.
편미분 방정식 (25)를 푸는 표준 방법은 푸리에 변환(Fourier transform)이다. 함수 $u(x,z)$의 푸리에 변환 $\widetilde{u}(x, p)$를 식 (25)에 넣어서 상미분 방정식으로 바꾼다.

                      (31a)

                      (31b)

치환 $U(x, z)$ = $e^{i k_0 x} u(x, z)$로 인해 $d^2/dx^2$은 거의 0이라서 식 (31b)의 해는 쉽게 얻어진다.

                      (32)

식 (32)를 역변환해서 임의 위치에 생기는 $u(x, z)$를 유도할 수 있다.

                  (33)

식 (33)은 $x, z$축을 분리한 형태라서, 전자파가 $x$축을 따라서는 프레넬 영역(Fresnel region)의 조건으로 전파되는 모습을 잘 보여준다. 그래서 식 (33)은 푸리에 분할 단계 알고리즘(Fourier split-step algorithm)을 적용한 전형적인 해법이다. 왜냐하면 $x, z$ 변수를 한꺼번에 풀지 않고, 문제를 분할해 처음 단계는 $z$축에 대한 푸리에 변환만 적용하고 그 다음 독립된 단계에서 $x$축에 대한 미분 방정식을 풀어 전체 문제를 해결하기 때문이다. 이러한 알고리즘을 포물형 방정식에 적용한 식 (33)은 간략히 분할 단계 포물형 방정식(split-step parabolic equation, SSPE)으로 불린다.
문제를 컴퓨터로 쉽게 풀기 위해 위치 $x$를 $x$ = $x_0 + s \Delta x$로 이산화해서 식 (33)이 $s$에 대한 재귀식(recurrence formula)이 되도록 한다.

                      (34)

여기서 $x_0$은 송신 안테나 위치, $s$ = $1, 2, \cdots$이다. 식 (34)는 재귀식이라서 $x$ = $x_0$부터 출발해 길이가 $\Delta x$만큼 증가할 때의 $u(x, z)$값을 재귀적으로 얻는다.

[참고문헌]
[1] J. R. Kuttler and G. D. Dockery, "Theoretical description of the parabolic approximation/Fourier split-step method of representing electromagnetic propagation in the troposphere," Radio Sci., vol. 26, no. 02, pp. 381–393, Mar.–Apr. 1991.
[2] 허준, 양준모, 박민규, 서영광, 박용배, "포물형 방정식을 이용한 레이다 표적 탐지의 굴절에 의한 오차 분석 연구," 한국전자파학회논문지, 제33권, 제5호, pp. 403–409, 2022년 5월.
[3] The Radio Refractive Index: Its Formula and Refractivity Data, Recommendation ITU-R P.453-14, Aug. 2019.
[4] 박명훈, 전우중, 김현승, 권세웅, 문현욱, 이기원, "대기 굴절률 변화에 따른 레이다 성능 변화에 관한 연구", 한국전자파학회논문지, 제32권, 제8호, pp. 743–750, 2021년 8월.
[5] G. Apaydin and L. Sevgi, Radio Wave Propagation and Parabolic Equation Modeling, Hoboken, NJ, USA: John Wiley & Sons, 2017.

댓글 없음 :

댓글 쓰기

욕설이나 스팸글은 삭제될 수 있습니다. [전파거북이]는 선플운동의 아름다운 인터넷을 지지합니다.