Chapter 08

위성 위치 계산

거리만 알아서는 위치를 못 구한다. “어디에 있는 위성까지의 거리”인지 알아야 한다. 위성은 자기 위치를 좌표로 보내지 않는다. 대신 궤도를 설명하는 값 16개(궤도력)를 보내고, 수신기가 원하는 순간의 위치를 직접 계산한다. GPS 규격은 그 계산법을 표 하나(표 20-IV)에 담았다. 이 장에서는 그 표를 한 줄씩 따라간다.

궤도력의 값들(√A, e, M₀, Δn, ω, Ω₀, Ω̇, i₀, IDOT, 보정 계수 6개, t_oe)이 무엇을 뜻하는지 안다.
규격 표 20-IV의 순서대로 위성의 지구 고정(ECEF) 좌표를 계산한다.
위성 시계 보정과 상대론 보정(Δt_r = F·e·√A·sin E_k)을 계산한다.
신호가 날아오는 70 ms 동안의 지구 자전을 보정하지 않으면 생기는 오차를 안다.

표 20-IV를 따라서

계산은 크게 다섯 단계다.

평균 운동과 평균 근점 이각: n = √(μ/A³) + Δn, M_k = M₀ + n·t_k. t_k는 궤도력 기준 시각 t_oe에서 지난 시간(주 경계를 넘으면 604,800초를 더하거나 뺌).
케플러 방정식: M_k = E_k − e·sin E_k를 반복으로 풀어 이심 근점 이각 E_k를 얻는다(규격: 최소 세 번 반복).
궤도면 안의 위치: 진근점 이각 ν_k, 위도 인수 Φ_k = ν_k + ω, 그리고 지구가 완전한 구가 아니어서 생기는 2차 조화 보정(C_us·C_uc, C_rs·C_rc, C_is·C_ic)으로 u_k, r_k, i_k.
승교점 경도: Ω_k = Ω₀ + (Ω̇ − Ω̇_e)·t_k − Ω̇_e·t_oe. 지구가 도는 만큼(Ω̇_e = 7.2921151467 × 10⁻⁵ rad/s) 빼서 지구 고정 좌표로.
회전: 궤도면 좌표 (r cos u, r sin u)를 i_k와 Ω_k로 돌려 x, y, z.

비행 계획서로 비행기 위치 맞히기

출발 시각, 항로, 속도가 적힌 비행 계획서가 있으면 레이더 없이도 지금 비행기가 어디쯤인지 계산할 수 있다. 바람에 밀린 만큼의 보정 항까지 있으면 더 정확하다. 궤도력이 위성의 비행 계획서다. 기준 시각(t_oe), 궤도 모양 6개 값, 그리고 보정 항들이다.

CALCULATOR

궤도력 → 위성 좌표, 한 줄씩

궤도력 기준 시각에서 지난 시간 t_k 이심률 e

그림: 위에서 본 지구(원)와 위성 궤도, 지금 위치(빨강). 궤도력 값은 실제 GPS 위성과 비슷한 크기로 만든 설명용 예시다(√A 5,153.65 m^½, M₀ 1.2, Δn 4.5 × 10⁻⁹, ω 0.82, Ω₀ −2.5, Ω̇ −8.1 × 10⁻⁹, i₀ 0.9615, IDOT 2 × 10⁻¹⁰(라디안 단위), C_uc −1.5 × 10⁻⁶, C_us 8 × 10⁻⁶, C_ic 6 × 10⁻⁸, C_is −1 × 10⁻⁷ rad, C_rc 230, C_rs −30 m, t_oe = t_oc 345,600 s, a_f0 1.2 × 10⁻⁴ s, a_f1 −3.4 × 10⁻¹² s/s). 해볼 것: ① t_k를 움직이면 위성이 약 12시간에 한 바퀴 돈다(그림은 지구가 도는 좌표라 궤적이 휜다). ② 이심률을 0.03으로 키우면 케플러 반복 횟수가 늘고, 상대론 보정이 커진다. 검증 스크립트는 같은 값을 RTKLIB의 eph2pos()에 넣은 결과와 mm 단위로 맞춰 본다.

계산식과 상수(μ = 3.986005 × 10¹⁴, Ω̇_e = 7.2921151467 × 10⁻⁵, π = 3.1415926535898, t_k의 주 경계 처리): IS-GPS-200N 표 20-IV(sheet 1~2, p.106~107) · 같은 계산의 구현: RTKLIB ephemeris.c eph2pos()

위성 시계와 상대론 보정

위성 시계도 GPS 시간과 조금씩 다르다. 관제국이 그 차이를 다항식(a_f0 + a_f1·(t − t_oc) + a_f2·(t − t_oc)²)으로 보내 준다. 여기에 3장에서 미리 맞춰 두지 못한 상대론 효과가 남는다. 궤도가 완전한 원이 아니어서(이심률 e), 위성이 지구에 가까울 때는 빠르고 낮은 곳에서 시계가 느리게, 멀 때는 반대로 간다. 이 주기적 효과가 Δt_r = F·e·√A·sin E_k, F = −4.442807633 × 10⁻¹⁰ s/m^½다. e = 0.01이면 최대 약 23 ns, 거리로 7 m다.

시계 보정 Δt_sv = a_f0 + a_f1(t − t_oc) + a_f2(t − t_oc)² + Δt_r, Δt_r = F e √A sin E_k, F = −2√μ/c² = −4.442807633 × 10⁻¹⁰: IS-GPS-200N 20.3.3.3.3.1(p.98) · RTKLIB은 같은 값을 −2√(μA)·e·sin E/c²로 계산(eph2pos() “relativity correction”)

날아오는 동안 지구가 돈다

위성 위치는 신호를 보낸 순간의 것이어야 하고, 좌표는 수신기가 받는 순간의 지구 고정 좌표여야 한다. 신호가 약 70 ms 날아오는 동안 지구는 Ω̇_e × 0.07 s ≈ 5.1 × 10⁻⁶ rad 돈다. 위성 좌표가 100~160 m 옮겨지는 셈이고, 그중 시선 방향 성분이 거리 오차가 된다. 동·서쪽 낮은 위성에서 약 30 m다. 그래서 보낸 순간의 위성 좌표를 그만큼 z축 둘레로 돌려 쓴다(사냐크 보정).

예측해 보기

신호가 위성에서 지상까지 약 70 ms 날아오는 동안, 적도의 수신기는 지구 자전으로 얼마나 움직일까?

적도의 자전 속도는 약 465 m/s.

465 m/s × 0.07 s ≈ 33 m. 그래서 송신 시각의 위성 좌표를 그 시간만큼 지구 자전 반대로 돌려 줘야 한다(사냐크 보정). 빼먹으면 위성 방향에 따라 수십 m가 틀린다. 아래 계산기에서 위성 방향을 바꿔 보자.

CALCULATOR

지구 자전 보정의 크기

위성 앙각 위성 방위

비행시간—

위성 좌표가 도는 거리—

거리 오차(보정 안 하면)—

수신기: 위도 37.5°, 경도 127°, 위성: 반지름 26,560 km. 해볼 것: 동쪽·서쪽 낮은 위성에서 거리 오차가 가장 크고(수십 m), 머리 위나 남북 방향에서는 작다. 계산: 보정 = Ω̇_e·(x_s·y_r − y_s·x_r)/c 를 거리로(RTKLIB geodist()와 같은 식).

핵심 정리

위성은 좌표가 아니라 궤도력(케플러 요소 + 보정 항)을 보내고, 수신기가 표 20-IV로 위치를 계산한다.
평균 운동 → 케플러 방정식(반복) → 진근점 이각 → 2차 조화 보정 → 지구 자전을 뺀 승교점 → 회전.
위성 시계 = a_f0 + a_f1Δt + a_f2Δt² + Δt_r, Δt_r = F·e·√A·sin E(e = 0.01이면 최대 약 7 m).
신호가 날아오는 70 ms의 지구 자전을 보정한다(사냐크). 안 하면 수십 m.

확인 퀴즈

케플러 방정식 M = E − e·sin E에서 수신기가 반복으로 구하는 것은?

M은 시간에 비례해 쉽게 구하고, E는 닫힌 식이 없어 반복한다.

Ω_k 계산에서 Ω̇_e·t를 빼는 이유는?

궤도면은 우주에 대해 거의 고정이고 지구가 그 아래서 돈다.

상대론 보정 Δt_r가 0인 궤도는?

Δt_r ∝ e. 원 궤도의 일정한 효과는 3장의 주파수 오프셋이 이미 처리했다.

사냐크(지구 자전) 보정이 필요한 이유는?

약 70 ms 동안 지구가 돌아 위성 좌표를 그만큼 돌려야 한다.