포물형 방정식(PE)에서 유효 굴절률의 새로운 근사 공식

Ⅰ. 서 론

장거리 전자파 전파를 분석하는 효율적인 수치 해석 기법으로 포물형 방정식(Parabolic equation)이 있다[1]-[3]. 여러 포물형 방정식 기법 중에서 푸리에 변환(Fourier transform)과 영역 이산화를 쓰는 분할 단계 포물형 방정식(SSPE, Split-Step Parabolic Equation)이 주로 활용된다[3]. 하지만 100 km 이상의 장거리 전송에는 파수를 설정하는 대기 굴절률(Atmospheric refractive index)이 큰 영향을 주기 때문에, 대기 조건의 측정값을 자세히 분석해서 대기 굴절률을 추정하고 고도별 대기 굴절률의 변화를 정확히 계산해야 한다. 레이다 운영 중에 대기 환경 변화로 인해 발생하는 성능 열화를 방지하기 위해 M&S(Modeling and Simulation) 방법을 이용하여 저고도 및 원거리 탐지 특성을 분석할 수 있다[4]. 혹은 미리 정보를 알고 있는 협조 표적의 정보와 라디오존데(Radiosonde)의 측정값을 사용하여 실제 대기 굴절률을 추정하는 기법도 제안되었다[5].

본 논문에서는 등각 사상에 속하는 뫼비우스 변환(Möbius transformation)을 활용해 곡률을 가진 지구 표면을 평평한 전파 모형으로 바꾼다[6]. 이때 발생하는 실제 고도와 변환된 고도의 오차 특성을 복소 평면에서 분석한다. 그후 맥스웰 방정식에 뫼비우스 변환과 자율 미분 방정식(Autonomous differential equation)을 적용함으로써, 편파별로 대기 굴절률 변화의 고차항까지 고려한 새로운 유효 굴절률(Effective refractive index) 공식을 제안한다[7]. 이 근사 공식을 바탕으로 라디오존데의 측정값을 비선형 최소 제곱법인 가우스–뉴턴 알고리즘(GNA, Gauss–Newton Algorithm)으로 추정한다[8]. GNA의 추정 함수로 설정한 대기 굴절률로 양방향 넓은 각 SSPE(2-way wide-angle Split-Step Parabolic Equation)를 계산해서 기존 고도별 대기 측정값 결과와 유사함을 보였다[3]. 본 논문의 방법론은 100 km 이상 범위에서 고도별로 복잡하게 변하는 대기상의 전자파 전파를 분석할 때 유용하다. 새로운 유효 굴절률 공식은 고도별 대기 굴절률의 변량까지 포함하기 때문에 복잡한 대기 변화로 인한 전자파의 굴절을 더욱 정확하게 예측할 수 있다. 또한 라디오존데의 많은 측정값을 쓰지 않고, GNA 기반으로 이 측정값을 2개의 변수로 구성한 추정 함수는 가변하는 대기 환경의 모형화에 효과적이다. 개선된 유효 굴절률과 GNA를 활용한 대기 굴절률의 예측은 장거리 전자파 통신 및 레이다 운영의 분석에 유용할 것이다. 특히 본 연구는 대기 환경이 장거리에 걸쳐 변하거나 수분이 많은 무선 통신 환경에서 레이다의 표적 추정 오차를 개선하거나 지하 공동구의 채널 특성을 분석할 때 쓸모가 많을 것이다[9][10].

Ⅱ. 뫼비우스 변환

맥스웰 방정식의 페이저(Phasor)를 위한 시간 약속(Time convention)은 e^-iωt로 가정한다. 그림 1에 보인 지구 표면상의 위치를 표현하기 위해 2차원의 극 좌표계(Polar coordinate system)를 도입한다.

Fig. 1.

Electromagnetic propagation over the Earth's surface represented in a spherical coordinate system

여기서 복소수 ξ,γ는 각각 y를 고정한 데카르트 좌표계 (x,z)와 ϕ가 일정한 구 좌표계 (r,θ)를 위한 복소 평면의 변수, R_e는 지구의 평균 반지름이며 6,371 km로 생각한다. 식 (1)에 가우스 광학(Gaussian optics)에 나오는 근축 근사(Paraxial approximation)를 적용하면, θ ≪ 1로 인해 cosθ ≈ 1, sinθ ≈ θ가 얻어진다. 여기에 더해 안테나의 설치 높이 h가 지표면 가까이 있다는 조건 h ≪ R_e에 따라 r ≈ R_e + h, x = r sinθ ≈ R_eθ로 근사화된다. 그러면 식 (1)은 θ,h로 다시 표현된다.

여기서 Δr = r - R_e ≈ h이다. 요건 θ ≪ 1 및 h ≪ R_e를 식 (2)에 대입하면, ξ ≈ 0이 되기 때문에 식 (1)에 의해 γ ≈ iR_e가 성립한다. 다음 단계로 휘어진 구 표면을 평평한 평면으로 만들기 위해 등각 사상(Conformal mapping)의 일종인 뫼비우스 변환(Möbius transformation)을 식 (3)에 도입한다[1][6].

여기서 각 계수는 ad - bc ≠ 0을 만족해야 한다. 식 (1)의 우변을 뫼비우스 변환의 형태로 만들려고 다음과 같은 근사를 채택한다.

여기서 식 (2)의 성과인 γ ≈ iR_e를 가정한다. 식 (4)는 식 (1)과 약간 다르므로, 식 (5)에 보인 참값에서 벗어나는 고도 오차 Δz를 복소 평면에서 정의한다.

여기서 z = h로 둔다. 고도 오차 Δz는 그림 1과 같은 둥근 지구 표면상 좌표 z를 뫼비우스 변환에 따라 평평한 표면의 좌표 z′으로 바꿀 때 발생하는 오차를 정량화한다. 그림 2는 뫼비우스 변환 근사로 인해 발생하는 거리에 대한 고도 오차 Δz의 변화를 보여준다. 안테나의 설치 높이 h에 따라 Δz가 비례하여 커진다. 또한 θ = 0°에서 Δz의 최대값이 식 (6)처럼 얻어진다.

Fig. 2.

Variation of height error Δz versus range x induced by the Möbius transformation

이 고도 오차 Δz는 h의 제곱에 비례하고 지구 반지름에는 반비례한다.

Ⅲ. 등각 사상과 유효 굴절률

구 좌표계로 표현한 맥스웰 방정식에 r ≫ 1 및 전방향 복사(Omnidirectional radiation) 조건 ∂/∂ϕ = 0을 부여하고 급속 하강 방법(Method of steepest descent)을 적용해서 수직 편파(VP, Vertical Polarization)가 만드는 자기장 H_ϕ(r,θ)의 표현식을 정의한다[1][3].

여기서 $$ ϵ\left(r,\theta \right)={ϵ}_0{ϵ}_r\left(r,\theta \right)$$, $$ k_0=\omega \sqrt{{\mu }_0{ϵ}_0}$$; $$ n\left(r,\theta \right)=\sqrt{{ϵ}_r\left(r,\theta \right)}$$는 대기 굴절률(Atmospheric refractive index)이다. 식 (7)로 나타낸 자기장 H_ϕ(r,θ)는 다음 파동 방정식을 만족한다. 이 파동 방정식은 맥스웰 방정식의 발산(Divergence)과 회전(Curl) 연산자를 구 좌표계에서 전개해서 구한다.

여기서 ∂_r = ∂/∂r, ∂_θ = ∂/∂θ이다. 식 (7)을 식 (8)에 대입해 ϕ축에 독립인 정규화 자기장(Normalized H-field) U(r,θ)의 2차원 파동 방정식을 만든다.

여기서 $$ {\nabla }_{r\theta }^2$$는 2D 라플라시안(Laplacian)이며 일반화 파수(Generalized wavenumber) Λ_rθ를 규정한다:

식 (10)에 나온 2계 편미분 기호는 $$ {\partial }_r^2={\partial }^2/\partial r^2$$, $$ {\partial }_{\theta }^2={\partial }^2/\partial {\theta }^2$$를 뜻한다. 식 (9)에 등각 사상인 식 (4)를 대입해서 x,z에 대해 정리한다.

여기서 좌표 변환에 필요한 미분은 식 (12)와 같다:

다시 식 (11)에 근축 근사인 $$ U\left(x,z\right)=e^{i{k}_{0}x}u\left(x,z\right)$$를 적용해서 포물형 편미분 방정식을 식 (13)처럼 유도한다.

여기서 $$ n_{\text{eff}}^2=1+n_d^2$$; u(x,z)는 축퇴 근축장(Degenerate paraxial field), n_eff는 VP에 대한 유효 굴절률(Effective refractive index), 유효 굴절률과 자유 공간 굴절률의 등가적 변화량 n_d는 아래처럼 계산된다:

등가적 변화량 n_d는 식 (4)의 뫼비우스 변환에 의한 파동 방정식의 공간 변환 정도를 가상적인 굴절률 차이로 치환한 결과이다. 식 (10)을 x,z 관점에서 좌표 변환하고 식 (14)에 대입함으로써 최종적으로 VP에 대한 n_eff의 근사 공식을 완성한다.

여기서 편미분을 간략화하기 위해 ∂_r ≈ ∂_z, ∂_θ ≈ r∂_x로 근사한다. 원역장 조건을 써서 x ≫ λ₀로 가정하면 식 (15)는 더욱 단순해진다.

여기서 식 (17)에 나온 Δn²은 Δn의 제곱이다. 식 (16)에 대해 z ≪ R_e, Δn ≈ 0을 선택하면, 기존에 사용하는 유효 굴절률 공식이 식 (18)처럼 유도된다[2].

수평 편파(HP, Horizontal Polarization)에 대한 유효 굴절률 m_eff도 비슷하게 유도할 수 있다. 먼저 HP가 생성하는 전기장 E_ϕ(r,θ)를 식 (19)와 같이 선택한다[3].

정규화 전기장(Normalized E-field) V(r,θ)는 식 (20)에 기술한 2차원 파동 방정식을 만족한다.

여기서 HP의 일반화 파수 Ξ_rθ를 식 (21)처럼 설정한다:

모든 결과를 종합하여 HP에 대한 유효 굴절률 m_eff를 식 (22)에 정확히 유도한다.

편파 VP와 HP에 대한 유효 굴절률 n_eff와 m_eff는 그림 1과 같은 지구 표면상의 전파 특성을 데카르트 좌표계로 편평화할 때 출현하는 파동 방정식의 왜곡을 등가적인 파동의 굴절률로 고친 양이다. 식 (16)과 동일한 근사 조건으로 m_eff를 더욱 간단히 표현한다.

VP로 인해 파생되는 Δn²의 성질을 파악하기 위해, Δn² = 0을 만족하는 해 n(z)를 자율 미분 방정식(Autonomous differential equation)으로 유도한다[7]. 자율 미분 방정식은 미분 계수 dn(z)/dz가 고도 z에는 독립이면서 대기 굴절률 n(z)에만 비례하는 조건을 가진다.

여기서 n₀ = n(z₀), a₀는 대기 굴절률의 상대 기울기(Relative slope of refractive index)이다. 대기 굴절률이 식 (24)처럼 변화하면 n_eff = m_eff이 되어서 식 (23)을 모든 편파용 유효 굴절률로 쓸 수 있다. 이어서 대기 굴절률의 고도 지수(Height exponent)인 q를 추가해 식 (24)를 일반화 대기 굴절률로 바꾼다.

식 (17)에 나오는 대기 굴절률의 구배(Gradient)를 식 (26)처럼 얻기 위해, 식 (25)를 z에 대해 미분한다.

대기 굴절률은 고도가 높아짐에 따라 항상 낮아지므로, ∂_zn은 음수가 되어야 해서 항상 qa₀ < 0을 만족해야 한다. 따라서 VP를 위한 유효 굴절률에 쓰이는 Δn이 식 (27)과 같이 간략히 결정된다.

대기 굴절률 n(z)는 1에 가까워서 n = 1에서 이격하는 정도를 명확히 나타내기 위해, 대기 굴절률의 N단위(N-unit)와 유효 굴절률의 M단위(M-unit)를 정의해 각각 사용한다. 식 (28)에 나오는 굴절률의 정의 N,M은 매우 작은 대기 굴절률 변동을 10⁶만큼 과장해서 표현하기 위해 도입한다.

특히 VP에서는 n(z)가 기존보다 작아지기 때문에 N단위 차이(N-unit difference) ΔN을 식 (29)처럼 활용한다.

N단위 차이 ΔN은 대기 굴절률 n을 기준으로 굴절률의 고도 변화율 ∂_zn을 반영한 Δn에 의한 굴절률의 변동량을 N단위로 나타낸다. 그림 3은 VP의 전자파 전파에 대해 a₀,q에 따라 ΔN이 바뀌는 성질을 제시한다. 대기 굴절률의 기울기 |qa₀|이 증가할수록 ΔN은 점점 증가해서 기존 n(z) 특성에서 멀어진다. 하지만 |a₀|이 10^-6 m^-1인 경우는 고도 지수 q가 달라져도 ΔN은 거의 0에 가깝다.

Fig. 3.

Dependence of N-unit difference ΔN on height z for vertically polarized waves (a) |a0| = 10-4 m-1 (b) |a0| = 10-6 m-1

운영 주파수 f를 조정할 때 발생하는 ΔN의 변동성은 그림 4에 나타낸다. HF(High Frequency) 대역에서는 ΔN이 z에 따라 유의미하게 증가하지만, HF보다 높은 주파수 대역에서는 ΔN이 0에 가까워서 Δn²을 무시할 수 있다.

Fig. 4.

Variability of N-unit difference ΔN with respect to operating frequency f

Ⅳ. 가우스–뉴턴 알고리즘

라디오존데(Radiosonde)로 매일 측정하는 대기 조건에는 고도, 습도, 대기압 등에 측정 오차가 포함되기 때문에 적절한 특성 추정이 필수적이다[2][4]. 통상적으로 선형 최소 제곱법(Linear least squares)이 측정 오류를 줄이는 최적 추정기(Optimal estimator)이지만, q ≠ 1인 식 (25)는 비선형 함수이다. 이때는 최소 제곱 행렬과 반복 연산에 기반한 가우스–뉴턴 알고리즘(GNA, Gauss–Newton Algorithm)이 유용하다[8]. 라디오존데의 측정값으로 식 (25)를 구성하는 a₀,q를 적합하게 결정하기 위해 식 (30)과 같은 GNA를 반복해서 계산한다.

여기서 $$ l=\mathrm{0,1},2,\cdots ,N_a-1$$, ; 재귀하는 제l번 변수 벡터(Variable vector)는 $$ {\overline{\beta }}^{\left(l\right)}={\left[a_0^{\left(l\right)},q^{\left(l\right)}\right]}^T$$, $$ {\mathit{J}}_{r,\beta }$$는 $$ \overline{r},\overline{\beta }$$에 대한 야코비 행렬(Jacobi matrix), $$ \overline{r}$$은 반복을 통해 0으로 수렴시키는 잉여 벡터(Residual vector)이다. 잉여 벡터 $$ \overline{r}$$ 및 야코비 행렬 $$ {\mathit{J}}_{r,\beta }$$의 원소는 다음처럼 생성한다.

여기서 식 (31)에 나오는 n_i는 z_i에서 측정한 대기 굴절률, 식 (32)는 야코비 행렬의 원소를 표현한다.

그림 5는 라디오존데 측정값과 가우스–뉴턴 알고리즘으로 추정한 N단위의 변화를 보여준다.

Fig. 5.

Estimation of the atmospheric refractive index using the Gauss–Newton algorithm (a) WMO Station ID 47122: z0 = 52 m, N0 =378.074 (b) WMO Station ID 47186: z0 = 235 m, N0 = 345.411

여기서 N₀은 z₀에서 결정한 N단위이다. 경기도 오산시(WMO Station ID 47122)와 제주도(WMO Station ID 47186)에서 잰 라디오존데의 대기 조건은 각각 2023년 7월 21일 및 2024년 5월 11일에 측정되었다. 여기서 WMO Station ID는 세계기상기구(World Meteorological Organization)에서 정한 전지구적 기상 관측소 식별자이다. 라디오존데의 날짜별 관측 결과는 미국 와이오밍 대학교 대기과학과(Department of Atmospheric Science, University of Wyoming, USA)에서 제공한다[11]. 라디오존데를 띄운 지역과 날짜가 다르더라도 GNA는 식 (25)를 기준으로 측정값을 잘 추정한다. 그림 5(a)와 5(b)에서 GNA로 획득한 최종값은 각각 a₀ = 0.00099543, q = -0.00011841 및 a₀ = 0.00040950, q = -0.00015275이다. 그림 5(a)에 나타낸 라디오존데의 고도별 측정값과 GNA를 이용한 대기 굴절률의 추정 함수를 각각 적용해 계산한 양방향 넓은 각 SSPE(2-way wide-angle Split-Step Parabolic Equation)의 결과를 그림 6(a)와 6(b)에 소개한다[3]. 양방향 SSPE 계산은 총 3회 시행되었다. 송신 안테나는 (x,z) = (0,600 m)에 위치하고, 안테나의 주빔 방향과 반전력 빔폭은 각각 -1° 및 1°이다. 산란파를 생성하는 지표면과 지형 구조는 [12]와 같이 PEC이며, 대기 굴절률을 제외한 나머지 계산 조건도 [12]와 동일하게 맞춘다. 다만 전방파(Forward wave)만 고려한 단방향(1-way) SSPE를 쓰는 [12]와 다르게 그림 6은 후방파(Backward wave)의 산란까지 포함된 양방향(2-way) SSPE로 계산되었다. 양방향 SSPE는 단방향 SSPE보다 더 정확한 근축장 분포를 생성한다. 모든 거리 x와 고도 z에 대해, 그림 6의 크기 분포는 거의 유사한 특성을 보이므로, GNA의 대기 굴절률 추정은 양방향 SSPE 처리에 유용하다.

Fig. 6.

Magnitude distributions of the degenerate paraxial field u(x,z) computed using the 2-way wide-angle SSPE for vertical polarization (a) Radiosonde measurement (b) Gauss–Newton algorithm

Ⅴ. 결 론

등각 사상인 뫼비우스 변환을 사용하여 휘어진 지구 표면을 평평한 지면으로 바꿀 때 발생하는 고도 오차의 최대값을 계산하였다. 이어서 뫼비우스 변환에 따라 맥스웰 방정식을 수정하여 유도한 유효 대기 굴절률의 개선된 공식을 수직 및 수평 편파에 대해 새롭게 제안하였다. 또한 라디오존데의 측정값을 비선형 최소 제곱법인 가우스–뉴턴 알고리즘으로 추정함으로써 대기 굴절률을 편리하게 다루는 방식도 제시하였다. 라디오존데의 이산적 측정값과 가우스–뉴턴 알고리즘의 연속적 추정 함수를 이용해 계산한 양방향 넓은 각 SSPE의 결과는 서로 유사한 특성을 보인다. 개선된 대기 굴절률 관계식과 가우스–뉴턴 알고리즘을 활용하여 향후 레이다 시스템의 주파수별 표적 인식 특성을 분석할 계획이다.

References

• 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 Science, Vol. 26, No. 2, pp. 381-393, Mar.-Apr. 1991. [https://doi.org/10.1029/91RS00109]
• G. Apaydin and L. Sevgi, "Radio Wave Propagation and Parabolic Equation Modeling", Hoboken, NJ, USA: John Wiley & Sons, pp. 40-42, 2017. [https://doi.org/10.1002/9781119432166]
• Y. H. Cho and K. I. Kwon, "An SSPE technique with impedance boundary conditions for enhanced stability and accelerated computations", The Journal of Korean Institute of Information Technology, Vol. 22, No. 11, pp. 61-69, Nov. 2024. [https://doi.org/10.14801/jkiit.2024.22.11.61]
• M.-H. Park, W.-J. Jeon, H.-S. Kim, S.-W. Kwon, H.-W. Moon, and K.-W. Lee, "Radar performance variation according to the change of atmospheric refractive index", The Journal of Korean Institute of Electromagnetic Engineering and Science, Vol. 32, No. 8, pp 743-750, Aug. 2021. [https://doi.org/10.5515/KJKIEES.2021.32.8.743]
• M.-G. Kim, Y.-S. Chung, S.-J. Ahn, C.-I. Hong, and K. I. Kwon, "Beam path modeling and atmospheric radio refractivity estimation based on co-operative target", The Journal of Korean Institute of Communications and Information Sciences, Vol. 49, No. 5, pp. 695-701, May 2024. [https://doi.org/10.7840/kics.2024.49.5.695]
• P. K. Kythe, "Handbook of Conformal Mappings and Applications", New York, NY, USA: CRC Press, pp. 51-56, 2019. [https://doi.org/10.1201/9781315180236]
• E. Kreyszig, "Advanced Engineering Mathematics", 10th Ed., Jefferson City, MO, USA: John Wiley & Sons, pp. 33-34, 2011.
• J. Solomon, "Numerical Algorithms: : Methods for Computer Vision, Machine Learning, and Graphics", Boca Raton, FL, USA: CRC Press, pp. 227-229, 2015. [https://doi.org/10.1201/b18657]
• M. Kim, Y. Kim, S. Park, S. Cho, C. Hong, and H. Kim, "Error analysis of radar target measurement due to atmospheric refractivity", Journal of The Institute of Internet, Broadcasting and Communication, Vol. 24, No. 6, pp. 57-62, Dec. 2024. [https://doi.org/10.7236/JIIBC.2024.24.6.57]
• B.-J. Lee and W.-S. Jung, "Characteristic analysis of wireless channels to construct wireless network environment in underground utility tunnels", Journal of The Institute of Internet, Broadcasting and Communication, Vol. 24, No. 3, pp. 27-34, Jun. 2024. [https://doi.org/10.7236/JIIBC.2024.24.3.27]
• Upperair Air Data, Department of Atmospheric Science, University of Wyoming, https://weather.uwyo.edu/upperair/sounding.html, . [accessed: Feb. 25, 2025]
• S. Ahn, J. Kim, K. I. Kwon, Y. H. Cho, and Y.-S. Chung, "Transparent boundary conditions for finite element method-based parabolic equation", Journal of Electromagnetic Engineering and Science, Vol. 25, pp. 1-3, Oct. 2024. [https://doi.org/10.26866/jees.2025.4.l.25]