Design of Orbit Simulation Tool for Lunar Navigation Satellite System

Abstract

Lunar Navigation Satellite System refers to a constellation of satellite providing PNT services on the moon. LNSS consists of main satellite and navigation satellites. Navigation satellites orbiting around the moon and a main satellite moves the area between the moon and the L2 point. The navigation satellite performs the same role as the Earth's GNSS satellite, and the main satellite communicates with the Earth for time synchronization. Due to the effect of the non-uniform shape of the moon, it is necessary to focus on the influence of the lunar gravitational field when designing the orbit simulation for navigation satellite. Since the main satellite is farther away from the moon than the navigation satellite, both the earth's gravity and the moon's gravity must be considered simultaneously when designing the orbit simulation for main satellite. Therefore, the main satellite orbit simulation must be designed through the three-body problem between the Earth, the moon, and the main satellite. In this paper, the orbit simulation tool for main satellite and navigation satellite required for LNSS was designed. The orbit simulation considers the environment characteristics of the moon. As a result of comparing long-term data (180 days) with the commercial program GMAT, it was confirmed that there was an error of about 1 m.

1. INTRODUCTION

Lunar Navigation Satellite System (LNSS)란 달 주위에 위성을 띄워 자체적인 PNT 서비스를 제공하는 항법 시스템을 의미한다. 근래에 들어 달탐사에 대한 관심이 전세계적으로 높아지면서 기존의 항법 기술의 한계를 극복하기 위한 방안으로 LNSS에 관한 연구가 활발히 진행되고 있다. 달에는 항법과 관련된 인프라가 전혀 없기 때문에 사용자는 맵 매칭 기법을 이용하여 위치를 계산하거나 지구와의 통신에 의존하여 항법을 수행할 수밖에 없다. 이러한 기존의 달 항법 기술들은 크게 두가지로 분류될 수 있으며, 달의 지형 정보를 활용하는 Landmark 기반 항법과 지구와의 통신을 활용하는 방법으로 나눌 수 있다. 먼저 지형을 활용한 항법의 대표적인 방법으로 Terrain Relative Navigation (TRN) 기법이 존재한다. TRN 기법은 탐사선이 달에 착륙하는 동안 카메라를 통해 획득한 이미지 정보와 미리 알고 있는 지형 정보를 비교하여 항법을 수행한다 (Johnson & Montgomery 2008). 비행중인 탐사선 뿐만 아니라 지표면을 이동하는 rover에도 적용 가능한 Landmark 기반 항법 기술을 개발하기 위해 NASA에선 Lunar Orbiter Laser Altimeter 데이터를 활용하여 능선과 바위 등을 AI 와 접목시키는 연구를 수행하고 있다. 이러한 landmark 기반 항법의 경우 사전에 고 정밀 지형 정보 가지고 있어야 한다는 점과 고해상도 카메라 및 고속 프로세서 등 고가의 장비를 필요로 한다는 단점이 존재한다. 지구와의 통신을 활용하는 방식에는 Deep Space Network (DSN) 신호를 활용하는 방법과 Global Navigation Satellite System (GNSS) 신호를 사용하는 방법 등이 존재한다 (Lanyi et al. 2008, Bhamidipati et al. 2021). 그러나 DSN은 채널 수에 제한이 있고 delay를 가지고 항법을 수행하기 때문에 실시간 항법엔 적합하지 않다. 또한 GNSS 위성 신호는 지구에 의한 음영 지역을 통과하기 전 후, 단 기간 동안만 수집할 수 있고 지구로부터 멀리 떨어진 곳에서 신호를 수집하기 때문에 신호 품질이 떨어지게 되며 위성의 극악한 기하학적 배치 문제가 가중되어 정확도를 저하 시키는 요인으로 작용한다. 또한 두 기법 모두 달의 뒷면에선 신호를 수신하지 못하는 한계 때문에 탐사 가능지역이 제한된다는 단점이 존재하게 된다.

LNSS는 기존의 항법 기술들의 한계를 해결하고자 연구중인 새로운 달에서의 항법개념이 이다. Fig. 1과 같이 LNSS는 주 위성과 항법 위성으로 이루어진다. 주 위성은 항법 위성들을 달까지 운반하는 모선의 역할을 수행함과 동시에 지구와의 통신을 위한 통신위성의 역할을 동시에 수행한다. 항법 위성은 GNSS 위성과 동일하게 항법 메시지를 전달하는 역할을 수행한다. 지구 중궤도에 한기씩 올릴 수 있는 GNSS 위성과 달리 LNSS의 항법 위성은 위성들을 지구에서부터 달까지 운반해 가야하기 때문에 동시에 여러 기를 운반할 수 있는 초소형급 사이즈의 위성을 활용하는 컨셉을 가지고 있다. 주 위성과 항법 위성은 각각 다른 운용목적을 가지고 있기 때문에 궤도 설계 시에도 고려해야 하는 요인이 다르다.

Fig. 1. Schematic picture of lunar navigation satellite system.

주 위성은 최초에 항법 위성들을 달까지 운반하는 모선 역할을 수행한다. 이후 달 궤도에 안착한 뒤부터는 지구와의 통신 위성 역할을 담당한다. 지구와의 통신을 통해 시각 동기화를 수행하고 다른 항법 위성들과 Inter Satellite Link를 통해 전체 시각 동기화를 수행한다. 따라서 주 위성 궤도는 지구와 원활한 통신을 위해 지구의 가시성이 높은 궤도를 선택해야 한다. 이때 달에는 달과 라그랑주 점 사이를 주기를 가지고 움직이는 3차원 형태의 halo orbit이 존재한다. 라그랑주 점은 지구와 달에 의한 중력과 위성의 원심력이 상쇄되어 나타나는 중력 평형점을 의미한다. 이러한 3차원 형태의 halo orbit들 중 커브가 거의 없는 형태를 갖는 궤도를 특별히 Near-Rectilinear Halo Orbit (NRHO)라고 부른다. NRHO 궤도는 지구가시성이 매우 좋은 궤도로 알려져 있다. 본 연구에선 LNSS의 주 위성 궤도를 L2 라그랑주 점과 달 사이를 움직이는 NRHO 궤도로 설계한다. 이때 NRHO 궤도는 달 중심으로부터 지구 중력이 우세한 영역까지 멀어지게 되기 때문에 시뮬레이션 설계 시 지구-달-위성의 3체 중력을 고려해야 한다.

항법 위성의 궤도는 궤도 요소의 변화가 적은 안정적인 궤도를 선택하는 것이 중요하다. 항법 위성에 활용하는 초소형 위성의 특성상 궤도의 안정성이 높을수록 궤도 수정을 위한 연료 소모를 최소화할 수 있으며 방송 궤도력을 생성하기에도 적합하다. 이때 항법 위성은 주 위성과 달리 상대적으로 달 근처를 움직이기 때문에 지구 중력보다 달 중력이 우세한 영향을 끼친다. 달의 경우 지구보다 형태적 불균형성이 강하며 이로 인한 영향은 spherical harmonics로 정의된 달 중력 포텐셜을 통해 고려할 수 있다 (Nie & Gurfil 2018). 지구 중력의 영향 또한 외란(third body perturbation)의 형태로 고려해야 한다. 본 연구에서는 앞서 설명한 각 궤도의 운용 환경 특성을 고려하여 궤도 설계 시 필요한 시뮬레이션 환경을 개발하였다. 본 연구에서는 앞서 설명한 각 궤도의 운용 환경 특성을 고려하여 궤도 설계 시 필요한 시뮬레이션 환경을 개발하였으며 상용 GMAT 프로그램과 비교하여 검증하였다.

2. LUNAR COORDINATE SYSTEM

2.1 Earth-Moon Geometry

지구는 Fig. 2와 같이 황도면을 기준으로 자전축이 약 23.5도 기울어진 채 태양 주위를 공전한다. 달의 공전궤도면은 황도면과 약 5.14도 기울어져 있으며 달의 자전축은 공전궤도면과 약 6.68도 기울어져 있다. 이때 달은 자전주기와 공전주기는 일치하는 조석 고정된 형태로 지구 주위를 움직이며 주기는 항성시로 약 27.3일이 소요된다. 조석 고정이란 지구에서 항상 달의 한쪽면만 보이는 현상을 의미한다. 이때 달의 공전궤도는 완벽한 원 형태가 아닌 이심률 0.0549의 타원 궤도이며 달의 자전축이 약 6.68도 기울어져 있기 때문에 달 칭동(libration) 현상이 발생한다.

Fig. 2. Lunar orbit plane relative to ecliptic plane.

2.2 Moon Libration

달의 칭동은 경도 칭동과 위도 칭동으로 나누어 설명할 수 있다. 먼저 경도 칭동은 달의 공전주기와 자전주기가 같지만 타원 궤도를 돌기 때문에 나타나는 현상이다. 달의 공전 속도는 원월점(apolune)에서 빨라지고 근월점(perilune)에서 느려지지만 자전속도는 일정하게 유지된다. 따라서 원월점과 근월점 근처를 이동할 때 상대적인 수평방향 움직임을 보인다. 위도 칭동의 경우 달의 자전축이 공전궤도면과 6.68도 기울어져 있기 때문에 발생하는 현상이다. 달의 자전축이 기울어진 채 지구를 공전하기 때문에 한 주기동안 남극 점과 북극 점 영역 중 많이 보이는 곳이 번갈아 나타나게 된다. 이로 인해 수직 방향의 상대적 움직임이 발생한다. 결과적으로 경도 칭동과 위도 칭동이 합해지면서 지구 입장에서 달을 바라볼 때 세차 운동 하는 듯한 상대적 움직임을 보인다.

2.3 MCI Coordinate Frame

궤도 설계를 위해선 달에서의 관성 좌표계를 정의해야 한다. 달의 관성 좌표계(Moon-Centered Inertial, MCI)는 지구의 관성 좌표계(Earth-Centered Inertial)인 International Celestial Reference System (ICRS) 좌표계와 유사하게 정의된다. MCI 좌표계는 달 중심으로부터 정의되며 z축 방향은 지구 시(Terrestrial time) 기준 2000년 1월 1일 12시(J2000)에 지구 자전축 방향으로 정의한다. x축은 같은 시점의 춘분점(vernal equinox) 방향으로 정의되며 y축은 z축과 x축의 외적 방향으로 결정된다.

2.4 MCMF Coordinate Frame

달 중심 고정 좌표계(Moon-Centered Moon-Fixed, MCMF)는 크게 두가지로 정의된다. Fig. 3과 같이 일반적으로 달의 위도 경도를 나타낼 때 이용하는 월면 좌표계는 달의 평균 자전축을 z축으로 정의하며 적도와 본초 자오선의 교점 방향을 x축으로 정의한다. 이때 앞서 언급한바와 같이 달은 칭동 현상으로 인해 지구를 바라보는 방향이 세차 운동하듯 변화하기 때문에 달의 본초 자오선 방향은 평균적으로 지구를 바라보는 방향으로 정의된다. 월면 좌표계는 좌표계를 정의하는 방법에 빗대어 Mean Earth (ME) 좌표계로 불리운다 (GSFC 2008). 또다른 MCMF 좌표계는 Principal Axis (PA) 좌표계이다. PA 좌표계는 달의 관성 모멘트로부터 정의된 주축(principal axis)으로 정의된다 (Standish & Williams 2012). ME 좌표계와 PA 좌표계는 일치하지 않기 때문에 달 표면에서 약 1 km까지 차이가 발생한다. NASA의 Jet Propulsion Laboratory은 행성들의 천체 궤도력을 제공한다. 이를 통해 좌표계간 변환할 수 있는 행렬을 계산할 수 있으며 본 연구에서는 여러 천체 궤도력 버전 중 DE421을 활용하였다. DE421에서 달 궤도는 Lunar Laser-ranging 데이터 피팅을 통해 서브 미터급의 정확도를 제공한다 (Folkner et al. 2009). ME와 PA 좌표계 사이의 변환은 식 (1)과 같이 회전 변환 행렬(C)을 통해 계산할 수 있다. ME와 PA 사이 각은 Table 1과 같이 DE421에서 고정된 상수로 제공한다 (Ye et al. 2018).

Fig. 3. Selenographic coordinate system.

Table 1. Transition angle between two MCMF frame (ME and PA).

	\(C_1\)	\(C_2\)	\(C_3\)
DE421	67.92″	78.56″	0.30″

 

\(\mathbf{C}=R_z (C_1 ) R_y (C_2 ) R_x (C_3 )\)                                                                                                                                              
\(=\left[ \begin{array}\cos ⁡C_1&\sin ⁡C_1&0\\ -\sin ⁡C_1&\cos ⁡C_1&0\\ 0&0&1\end{array}\right]\left[ \begin{array}\cos ⁡C_2&0&-\sin ⁡C_2\\ 0&1&0\\ \sin ⁡C_2&0&\cos ⁡C_2 \end{array}\right]\left[ \begin{array}1&0&0\\ 0&\cos ⁡C_3&\sin ⁡C_3\\ 0&-\sin ⁡C_3&\cos ⁡C_3 \end{array}\right]\)                              
\(r ⃗_{PA}=\mathbf{C}\cdot r ⃗_{ME}\)                                                                                                                                                                   (1)

특정 시점의 MCI와 PA 좌표계 사이의 오일러각  \(ϕ_m,θ_m,ψ_m\)은 DE421 천체 궤도력을 통해 계산할 수 있다. 이를 이용하여 식 (2)와 같이 MCI와 PA 좌표계 사이 변환 행렬(L)을 계산할 수 있다 (Taylor et al. 2010).

\(\mathbf{L}=R_z (-ϕ_m ) R_y (-θ_m ) R_x (ψ_m )\)                                                                                                                                        
\(=\left[ \begin{array}\cos ⁡ϕ_m&-\sin ⁡ϕ_m&0\\ \sin ⁡ϕ_m&\cos ⁡ϕ_m&0\\ 0&0&1\end{array}\right]\left[ \begin{array}\cos ⁡θ_m&0&\sin ⁡θ_m\\ 0&1&0\\ -\sin ⁡θ_m&0&\cos ⁡θ_m \end{array}\right]\left[ \begin{array}1&0&0\\ 0&\cos ⁡ψ_m&-\sin ⁡ψ_m\\ 0&\sin ⁡ψ_m&\cos \end{array}\right]\)                              
\(r ⃗_{MCI}=\mathbf{L}\cdot r ⃗_{PA}\)                                                                                                                                                                       (2)

3. ORBIT DYNAMICS

3.1 Main Satellite Orbit

주 위성의 후보 궤도인 NRHO 궤도는 달과 L2 점 사이에 존재하며 근월점에선 달의 영향이 우세하고 원월점에선 지구의 영향이 우세한 특징을 갖고 있다. 따라서 NRHO 궤도 시뮬레이션을 설계하기 위해선 지구-달-위성의 3체 운동방정식을 고려해야한다 (Zimovan et al. 2017). 위성의 질량은 지구와 달의 질량 대비 매우 작기 때문에 위성이 지구와 달에 미치는 영향은 무시 가능하다고 가정할 수 있다. 이때 지구와 달이 공통 질량 중심점을 기준으로 원 궤도를 움직인다고 가정하고 세운 3체 운동방적식을 원형제한삼체문제(Circular Restricted Three-Body Problem, CR3BP)라고 부른다 (Howell & Breakwell 1984). CR3BP 문제를 풀기위해선 지구와 달이 질량중심점(barycenter)을 중심으로 하는 새로운 회전 좌표계(barycentric rotating reference frame)를 설정해야 한다. Fig. 4는 지구와 달의 질량 중심점을 원점으로 하는 회전 좌표계를 의미한다. 이때 좌표계는 달의 공전 주기와 동일한 주기의 각속도를 가지고 z축으로 회전한다. 따라서 해당 좌표계에선 달과 지구가 항상 x축 위에 존재하게 된다. Fig. 5와 같이 지구와 달사이 거리를 1로 무차원화 시키면 질량중심점으로부터 지구와 달의 거리는 둘의 질량 비(\(\mu\)) 를 통해 나타낼 수 있다. 지구의 질량을 \(m_1\) 달의 질량을 \(m_2\)라 하였을 때 둘의 질량 비는 \(\mu = m_2/(m_1 + m_2)\)로 정의할 수 있다. 또한 달의 공전 각속도(rad/sec)를 이용하여 시간에 대해 무차원화 시키면 다음과 같은 3체 운동 방정식을 얻을 수 있다 (Vallado 2013).

Fig. 4. Position of the Earth and the Moon on the barycentric coordinate system.

Fig. 5. Definition of the CR3BP vector in the barycentric rotating reference frame.

\(x ̈=2y ̇+x-{1-\mu \over r_1^3} (x+\mu)-{\mu \over r_2^3} (x-1+\mu)\)                                                                 (3)

\(y ̈=-2x ̇+y-{(1-\mu )y\over r_1^3 }-{\mu y\over r_2^3 }\)                                                                                         (4)

\(z ̈=-{(1-\mu )z\over r_1^3}-{\mu z\over r_2^3}\)                                                                                                     (5)

where \(r ̄_1=[x+\mu ~~ y ~~ z]\), \(r ̄_2=[(x+\mu -1~~ y~~ z)]\)

여기서 \(r_1\)과 \(r_2\)는 각각 \(\bar{r}_1, \bar{r}_2\) 벡터의 norm 값을 의미한다. CR3BP에는 무수히 많은 주기 해(periodic solution)이 존재한다 (Grebow 2006). 그러나 비선형 방정식의 식 (3-5)는 closed form 형태의 해를 얻을 수 없다. 따라서 다양한 방식의 특수 해(particular solution)를 구하기 위한 연구가 진행되었고 본 연구에선 NRHO와 관련된 선행 연구의 초기 값을 활용하여 궤도 시뮬레이션에 사용하였다 (Zimovan 2017). 주 위성궤도의 위치와 속도를 \(\bar {\mathbf{x}}=[x~~y~~z~~\dot{x}~~\dot{y}~~\dot{z}]\)  벡터로 정의하게 되면 \(\dot{\bar {\mathbf{x}}}\)을 식 (3-5)을 통해 계산할 수 있다. 이를 8차 Runge-Kutta 기법을 통해 수치 적분하여 주 위성 궤도 전파를 수행하였다.

3.2 Navigation Satellite Orbit

항법 위성은 주 위성 대비 비교적 달 근처에 위치하기 때문에 달 중력의 영향을 가장 우세하게 받는다. 선행 연구에 따르면 항법 위성이 위치하게 고도에서는 달과 위성 사이 2체 중력이 가장 큰 영향을 미치며 그 외의 불 균일한 중력장에 의한 영향이나 지구 중력의 영향은 식 (6)과 같이 외란의 형태로 고려할 수 있다 (Nie & Gurfil 2018).

\(a ⃗=a ⃗_{Two~Body}+a ⃗_{Geopotential}+a ⃗_{Third~body}\)                                                                 (6)

달에서의 중력 포텐셜은 식 (7)과 같이 구면조화함수로 모델링할 수 있다. 

\(U(r,λ,ϕ)={μ \over r} \left[1+\sum\limits_{n=2}^∞\left({R\over r}\right)^n J_n P_{n0} (\sin ⁡ϕ ) +\sum\limits_{n=2}^∞ \sum\limits_{m=1}^n \left({R\over r}\right)^n \{C_{nm} \cos ⁡m λ+S_{nm} \sin ⁡m λ\} P_{nm} (\sin ⁡ϕ )\right]\)(7)

구면조화함수는 르장드르 다항식 \(P_{nm}\)과 계수 \(J_n, C_{nm}, S_{nm}\)으로 구성된다. 각각의 계수들은 달의 형태적 특성을 반영한다. 식 (7)에서 1은 구 형상을 의미하고 \(J_n\)은 회전 축에 대해 회전대칭 형상(zonal), \(C_{nm}\)은 체커 형상(tesseral), \(S_{nm}\)은 축 대칭 형상(sectoral)를 의미한다. Table 2는 NASA에서 제공하는 중력 포텐셜 모델 GRGM 900C 계수를 나타낸다 (Klokočník et al. 2022). 해당 모델은 NASA의 Goddard Space Flight Center (GSFC)의 중력 행성 임무인 Gravity Recovery and Interior Laboratory (GRAIL)임무를 통해 획득한 달 중력 모델로 2011년 9월에 발사한 두 대의 GRAIL 위성이 달 저궤도를 이동하며 획득한 데이터를 이용하여 제작되었다.

Table 2. Spherical harmonics coefficients of lunar gravity potential model.

GRGM 900C
\(J_n\)	\(n\) = 2 ⁝	-9.0886616361343905e-05 ⁝
\(C_{nm}\)	\(n\) = 2, \(m\) = 1 \(n\) = 2, \(m\) = 2 ⁝	-1.9495932361259600e-10 3.4673334085426700e-05 ⁝
\(S_{nm}\)	\(n\) = 2, \(m\) = 1 \(n\) = 2, \(m\) = 2 ⁝	1.0733325387331199e-09 1.0908462869835500e-10 ⁝

 

\(a ⃗_{Geopotential}=∇U(r,λ,ϕ)={∂U\over ∂r} e ̂_r+{1\over r}{∂U\over ∂ϕ} e ̂_ϕ+{1\over r \cos ⁡ϕ}{∂U\over ∂λ} e ̂_λ\)                                           (8)

식 (8)과 같이 중력 포텐셜의 그래디언트를 계산하면 달 중력장에 의한 가속도를 계산할 수 있다. 이때 식 (7)의 계수들이 모두 0인 경우 two-body에 의한 가속도( \(a ⃗_{Two~Body}\))를 의미하며 \(C_{nm}, S_{nm}\) 계수들에 해당하는 항은 달의 형태적 불 균일성으로 인한 외란(\(a ⃗_{Geopotential}\) )을 의미한다.

앞서 언급한바와 같이 항법 위성은 달 중력이 가장 우세하기 때문에 지구 중력의 영향은 주 위성에 비해 작지만 무시할 수 없는 수준으로 작용하기 때문에 지구 중력의 영향은 식 (9)와 같이 외란(\(a ⃗_{Third~body}\) )의 형태로 고려하게 된다.

\(a ⃗_{Third~body}=μ_{Earth}\left( {r ⃗_{Sat-Earth}\over r_{Sat-Earth}^3}-{r ⃗_{Moon-Earth}\over r_{Moon-Earth}^3}\right)\)                                                                 (9)

4. SIMULATION

주 위성과 항법 위성의 궤도 특성을 고려하여 궤도 전파를 수행하기 위한 각각의 시뮬레이션 툴을 설계하였다. 주 위성 궤도는 지구-달-위성의 3체 운동방정식을 고려하여 시뮬레이션 하였고 항법 위성 궤도는 달과 위성의 2체 운동방정식으로부터 불 균일한 달 중력장의 영향과 지구 중력을 외란으로 고려하여 시뮬레이션 하였다. 또한 항법 위성 궤도 시뮬레이션의 경우 상용 소프트웨어와 비교를 통한 검증을 수행하였다. 이때 상용 소프트웨어는 NASA에서 개발하고 오픈소스로 제공하는 GMAT 프로그램을 사용하였다.

주 위성의 궤도에 해당하는 NRHO 궤도는 CR3BP를 통해 무수히 많은 초기 조건 해를 얻을 수 있다. 이때 주 위성 궤도 전파시에는 NRHO 선행연구의 초기 조건을 이용하였으며 이는 Table 3과 같다 (Zimovan 2017). Figs. 6-8에 나타난 좌표계는 앞서 설명한 지구와 달의 질량 중심을 기준으로 달의 공전 속도와 동일한 각속도로 회전하는 회전 좌표계(barycentric rotating reference frame)를 의미한다. Fig. 6을 통해 NRHO 궤도군이 회전 좌표계 상에서 달과 L2 라그랑주 점 사이에 존재하는 것을 확인할 수 있다. Fig. 7은 y 방향에서 바라본 NRHO 궤도를 나타낸다. 달에 가까울수록 경사각이 높고 장반경은 달에서 L2 라그랑주 점으로 갈수록 점점 커지다가 줄어드는 것을 확인할 수 있다. Fig. 8은 x 방향에서 바라본 NRHO 궤도를 나타낸다. 이때 회전 좌표계는 달의 공전속도와 동일하게 회전하기 때문에 x축 방향에서 보았을 때 항상 달과 지구가 원점에 위치하게 된다. 이를 통해 NRHO 궤도에선 지구가시성이 높다는 것을 확인할 수 있다.

Table 3. Initial position and velocity of near rectilinear halo orbit.

\(x_0\) (km)	\(y_0\)(km)	\(z_0\)(km)	\(\dot{x}_0\)(km)	\(\dot{y}_0\)(km)	\(\dot{z}_0\)(km)
391808.96404 395083.47890 398340.34903 401532.01031 404650.40995 407705.32132 410964.34563 413951.41424 416948.17243 419991.86431 423131.57970 426434.09349 430302.36335 434269.91566 438741.26953 443687.90017 448581.84841 452299.85900 453933.39986	0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0	-69242.90379 -71423.21828 -73205.80678 -74663.49348 -75837.44224 -76744.69466 -77423.57581 -77747.17027 -77737.02288 -77324.63779 -76403.71956 -74810.22741 -72024.51252 -68031.50947 -62003.58623 -53037.39549 -40355.61358 -23889.07424 -1258.38952	0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0	-0.0948647458 -0.1127233517 -0.1290985769 -0.1438665551 -0.1571073553 -0.1689557414 -0.1803686386 -0.1896956286 -0.1979245392 -0.2050691756 -0.2110842836 -0.2158359883 -0.2192204355 -0.2200845466 -0.2175763789 -0.2095975208 -0.1936924761 -0.1706478547 -0.1521837491	0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

 

Fig. 6. NRHO group shown on the barycentric rotating reference frame.

Fig. 7. NRHO group viewed from the y direction on a barycentric rotating reference frame.

Fig. 8. NRHO group viewed form the x direction on a barycentric rotating reference frame.

항법위성 궤도의 궤도 전파 시에는 달의 중력장 모델과 지구 중력을 외란으로 고려하여 시뮬레이션을 설계하였다. 달 중력장 모델의 경우 GRGM 900C 모델을 이용하였다. 해당 모델은 상용 프로그램인 GMAT에서 사용하고 있는 중력 포텐셜 모델이다. 본 연구에선 GMAT과 비교 검증을 위해 동일한 버전의 중력 포텐셜 모델을 활용하였다. 이때 항법 위성 궤도의 후보군중 하나인 ELFO 궤도를 생성하여 GMAT과의 비교를 수행하였다. 일반적으로 ELFO 궤도를 설계하는 1,500 km 이상의 고고도 영역에서는 달의 2체 중력 다음으로 지구중력의 영향이 우세하며 그 다음으로 달 중력 포텐셜의 2차항에 해당하는 J2와 C22에 의한 섭동력이 우세한 영향을 갖는다 (Nie & Gurfil 2018). 따라서 본 연구에서는 ELFO 궤도 시뮬레이션 시 중력장 포텐셜에 의한 섭동력 중 우세한 영향을 끼치는 2차항까지 고려하여 GMAT과 비교하였다. 또한 3체 운동 방정식을 직접 고려하는 주 위성 궤도 시뮬레이션과 달리 항법 위성 궤도에서는 달 중력 다음으로 지구 중력이 우세하기 때문에 지구 중력은 외란으로 고려해 주었다. 지구 중력 고려시에는 천체 궤도력 DE421을 활용하였다. 시뮬레이션은 두가지 경우에 대해 분석을 진행하였다. 첫번째 경우는 Table 4의 ELFO 궤도 선행연구의 초기조건을 활용하여 한 주기 동안의 결과를 비교하였다 (Iiyama et al. 2022). 이때 선행 연구의 조건의 경우 궤도 수정 없이 단순 장기간 전파 시 달에 충돌하는 문제가 발생하기 때문에 추가적인 경사궤도에 대해 장기간 전파한 결과를 비교하였다.

Table 4. Initial Keplerian orbit elements of elliptical frozen orbit.

\(a\) (km)	\(e\)	\(i\) (°)	\(\Omega\) (°)	\(\omega\) (°)	\(M_0\) (°)
5053.7	0.6067	63.4	0.0	270.0	0.0

 

Fig. 9는 본 연구에서 개발한 항법 위성 시뮬레이션과 GMAT을 이용하여 선행 연구의 초기조건으로부터 각각 한 주기 동안 전파시킨 ELFO 궤적을 MCMF 좌표계상에 나타낸 결과이다. 이때 Fig. 10에서 각 축 별로 위치 오차를 살펴보게 되면 dm 수준의 차이를 갖는 것을 확인할 수 있다. Fig. 11은 각 궤적의 케플러 궤도 요소를 나타낸 결과이며 두 궤적의 궤도요소가 서로 동일한 개형을 갖는 것을 확인할 수 있다. 이때 궤도요소의 차이를 나타낸 결과는 Fig. 12와 같다. 각 궤도 요소별 오차 수준을 살펴보면 Fig. 11의 결과 값 대비 10^-6~10^-8 수준으로 무시 가능한 수준임을 확인할 수 있다.

Fig. 9. ELFO trajectory in MCMF coordinate system.

Fig. 10. Position error in MCMF coordinate system.

Fig. 11. Keplerian elements of navigation satellite orbit simulation and GMAT.

Fig. 12. Difference in Keplerian elements of navigation satellite orbit simulation and GMAT.

장기간 데이터 결과 비교를 위해 Table 5와 같이 새로운 궤적을 생성하여 90일간 궤도 전파를 수행하였다. Fig. 13은 생성된 궤적의 GMAT과의 위치 오차를 나타낸 그래프다. 이때 위치 오차가 장주기와 단주기가 결합되어 진동하는 형태를 가지며 각 축 별 오차 수준이 m 수준을 갖는 것을 확인할 수 있다. Fig. 14는 동일 궤적의 궤도 요소에 대해 GMAT과 차이를 나타낸 그래프이며 여전히 무시 가능한 수준임을 확인 할 수 있다. 이를 통해 본 연구에서 개발한 툴을 이용하여 초기 항법 위성 궤도 설계에 활용할 수 있을 것으로 판단된다.

Table 5. Initial Keplerian orbit elements for long term data comparison.

\(a\) (km)	\(e\)	\(i\) (°)	\(\Omega\) (°)	\(\omega\) (°)	\(M_0\) (°)
8000	0.6	63.4	0.0	90.0	0.0

 

Fig. 13. Position error in MCMF coordinate system for long term data (90 days).

Fig. 14. Difference in Keplerian elements of navigation satellite orbit simulation and GMAT for long term data (90 days).

5. CONCLUSION

본 연구에선 LNSS의 궤도 설계를 위한 시뮬레이션 툴의 제작 및 검증을 진행하였다. LNSS의 경우 지구와 통신 위성 역할을 수행하는 주 위성과 항법 메시지를 전달하는 항법 위성들로 이루어져 있다. 먼저 주 위성의 경우 지구와의 통신 위성 역할을 담당하기 때문에 지구와 통신이 원활할 수 있도록 지구 가시성이 높은 궤도를 설계해야 한다. 이때 여러 달 궤도들 중에는 지구 가시성이 높은 NRHO 궤도가 존재한다. NRHO 궤도는 원월점 부근에서 달 중력보다 지구 중력이 우세한 영역까지 멀어 지기 때문에 지구-달-위성의 3체 운동 방정식을 고려해야 한다. 따라서 주 위성 궤도 시뮬레이션 설계 시 3체 운동 방정식을 통해 툴을 제작하였다. 이를 통해 선행 연구의 NRHO 초기 조건을 활용하여 궤도 전파를 수행하였을 때 주기성을 가지면서 지구 가시성이 높은 궤도를 획득할 수 있는 것을 확인하였다. 항법 위성은 달의 중력이 우세한 영역을 움직이는 궤도를 이용한다. 따라서 항법 위성 궤도 시뮬레이션 설계 시 달과 위성의 2체 방정식에 달의 불규칙한 형상으로 인한 중력 포텐셜의 영향과 지구 중력의 영향을 외란으로 고려하였다. 달의 중력 포텐셜의 경우 NASA에서 제공하는 GRGM 900C 모델을 이용하여 2차항까지 고려하였으며 지구 중력은 DE421 천체 궤도력을 이용하였다. 이를 상용 프로그램인 GMAT의 결과와 비교를 통해 설계한 시뮬레이션 툴의 검증을 진행하였다. 검증에 사용한 항법 위성 궤도의 초기 값은 먼저 선행연구의 ELFO 초기 조건을 활용하였다. 한 주기의 궤도 전파 시뮬레이션을 통해 설계한 궤도 시뮬레이션과 상용프로그램 GMAT이 유사한 결과를 얻는 것을 확인하였다. 이후 고도와 이심률을 변화시킨 새로운 타원 경사 궤도를 통해 90일의 장기간 데이터를 비교한결과 GMAT과의 오차가 m 수준으로 발생하는 것을 확인하였다. 이를 통해 추후 LNSS 궤도 설계 시 본 연구에서 개발한 궤도 전파 시뮬레이션 툴을 활용할 수 있을 것으로 기대된다.

AUTHOR CONTRIBUTIONS

Conceptualization, H.J., methodology, H.J. and J.P., software, J.P., H.J and J.S., validation, H.J. and J.P., formal analysis, H.J. and J.P., investigation, J.P., H.J. and J.S., resources, H.J., and J.P., data curation, H.J., J.S., and M.K. writing—original draft preparation, H.S., and C.K., writing—review and editing, H.J.; visualization, H.J., supervision, C.K.

CONFLICTS OF INTEREST

The authors declare no conflict of interest.