Kalman Filter-Based Ensemble Timescale with 3- Hydrogen Masers

Abstract

A Kalman filter algorithm is used for the generation of an ensemble timescale with three hydrogen masers maintained in KRISS. Allan deviation curves of three pairs of clocks were obtained by a three-cornered hat method and were used as reference curves for determination of parameters of the Kalman filter-based timescale. The ensemble timescale equation of a 3-clock system was established, and the clocks' phases estimated by the Kalman filter were used as the prediction time of each clock in the equation. The weight of each clock was determined inversely proportional to the Allan variance calculated with the clocks' phases. The Allan deviation of the weighted mean was 1.2×10^-16 at the averaging time of 57,600 s. However when we made fine adjustments of the clocks' weight, the minimum Allan deviation of 2×10^-17 was obtained. To find out the reason of the great improvement in the frequency stability, additional researches are in progress theoretically and experimentally.

1. INTRODUCTION

Global Navigation Satellite System (GNSS)를 이용하여 정밀한 위치 추정과 항법이 가능한 것은 현대에 들어 정확한 시간을 생성하고 공급하는 것이 가능해졌기 때문이다. 정확한 시간 생성을 위해 제일 먼저 필요한 것은 정확하고 안정된 원자시계이다. 1950년대에 개발이 시작된 세슘원자시계를 비롯하여, GNSS 위성에 많이 탑재되는 루비듐원자시계, 단기 안정도에서 탁월한 성능을 가지는 능동형 수소메이저 등, 많은 종류의 원자시계들이 개발되어 상용화되었다 (Lee 2018). 하지만 이러한 원자시계들은 어느 한 순간에 고장이 나거나 수명이 다할 수 있다. 그러므로 위성항법시스템에서는 여러 대의 원자시계들을 모아서 앙상블(ensemble)을 구성하고, 하나 또는 그 이상의 시계가 정상적인 동작을 하지 못하더라도 정확하고 안정도와 신뢰도가 높은 시간눈금(timescale)을 만드는 것이 요구된다. 이런 앙상블 시간눈금을 생성하려면 시계들 사이에서 측정된 시간차이 데이터와 이론적인 시계모델을 바탕으로 시계들을 서로 결합하는 계산 알고리즘이 필요하다. Global Positioning System (GPS)의 경우, Brown(1991)은 “GPS composite clock”이라는 이름의 앙상블 시간눈금생성 이론을 확립했으며, 현재 GPS의 시스템 시간 생성에 적용되고 있다. 이 이론은 위성에 탑재된 시계로 측정한 GPS 항법신호의 송신 시각과 지상에 있는 시계로 측정한 수신 시각 사이의 시간차이 측정 데이터를 이용한다. 이 측정 데이터를 시계 모델과 결합하여 개별 위성에 탑재된 시계와 composite clock, 즉 앙상블 시간눈금과의 시간차이를 추정하는 데 칼만(Kalman) 필터 기반의 알고리즘이 사용된다 (Galleani & Tavella 2010).

칼만 필터는 1960년대 초에 R.E. Kalman이 처음으로 개발한 추정 기법으로, 칼만필터를 앙상블 시간눈금 생성에 처음으로 응용한 것은 1987년에 Jones와 Tryon이었고, Greenhall (2013)은 Jones-Tryon 모델을 발전시켜 몇가지 변형된 칼만필터 기반의 시간눈금 생성 알고리즘을 개발했다. 그런데 이들은 위성이 아니라 실험실에서 관리되고 있는 시계들로부터 시간눈금을 생성하는 것을 목표로 했다 (Greenhall 2006). 이 외에도 KAS-2는 미국 특허를 받은, 좀더 정교하게 다듬은 칼만 필터 기반의 알고리즘이다(Stein 1992). Levine은 이 알고리즘들을 비교 분석하여 각각의 장단점과 특성을 발표했다(Levine 2012). 그리고 최근에는 칼만필터를 사용하여 UTC(세계협정시)를 생성하는 새로운 방법이 소개되었다 (Parisi & Panfilo 2016, Panfilo & Arias 2019).

칼만필터를 통해 생성되는 앙상블 시간눈금은 개별 시계들이 갖는 주파수안정도보다 더 안정된 시간눈금을 얻을 수 있다는 장점이 있다. 이 보다 더 중요한 것은 어느 하나의 우수한 시계에 의존하는 것보다 덜 우수하더라도 여러 대의 원자시계를 사용함으로써 시간 생성의 연속성과 신뢰성을 유지할 수 있다는 것이다. 칼만필터는 안정된 시간신호를 만들어내지만 그 출력물이 나타내는 시각이 어떤 기준 시간과 일치하도록 조정(steering)되어야 한다.

예를 들면, GPS의 시스템 시간(이하, GPS 시간)은 미국 워싱턴 D.C.에 위치한 해군관측소(USNO)에서 생성하는 시간눈금인 UTC(USNO)와 시간이 일치하도록 조정된다(Kaplan & Hegarty 2006). 단, 윤초가 적용되지 않는 GPS 시간과 윤초가 적용되는 UTC와의 시간 차이는 위성에서 방송하는 항법신호 정보에 포함되어 전송되고, 이를 수신한 GPS수신기에서 보정한다. USNO에서는 UTC(USNO)를 생성하기에 앞서 USNO에서 관리하고 있는 여러 원자시계들을 이용하여 앙상블 시간눈금인 TA(USNO)를 먼저 생성한다. 이를 위해 BIPM에 등록되어 있는 USNO의 시계 수는 2020년 5월 현재 총 91대이다. 그 중 53대가 세슘원자시계이고, 34대가 수소메이저이며, 4대의 루비듐원자분수시계도 포함되어 있다 (BIPM Database 2020).

유럽우주국(ESA) 주관으로 2016년 12월부터 서비스 중에 있는 Galileo 시스템의 경우 그 시스템 시간(GST)은 유럽 여러 나라의 표준연구기관(독일 PTB, 이태리 IT, 프랑스 OP 등)에서 유지 및 관리하는 UTC 시간에 일치하도록 조정된다 (Lewandowski & Arias 2011). 적어도 3개 이상의 표준연구기관이 매일 Galileo Time Service Provider (GTSP)로서 정확한 시간을 제공하므로 이들과 매일 접속하는 지상의 위성통제센터에는 Precise Time Facility (PTF)로서 2+1 대의 능동형 수소메이저와 3+1대의 고성능 세슘원자시계를 갖출 것을 제안했다 (Bedrick et al. 2004).

그리고 중국에서 운영 중인 BeiDou 시스템 시간(BDT)의 경우는 중국 National Timing Service Center (NTSC)에서 유지하는 UTC(NTSC)에 일치하도록 조정된다 (Han et al. 2011). GNSS의 시스템 시간과 연계되어 있는 이들 표준연구기관(lab)에서 생성하는 UTC(lab)의 시간은 국제기구인 BIPM(국제도량형국)에서 생성하는 UTC와 수 나노초(ns) 이내에서 일치하도록 조정되고 있다. 결국 항법시스템을 운영하는 지역 또는 국가의 UTC(lab)을 정확하고 안정되게 생성하는 것이 GNSS 시간의 정확도를 높이는 출발점이다 (Lewandowski & Arias 2011).

본 논문에서는 한국표준과학연구원(KRISS)에서 유지·관리하고 있는 3대의 수소메이저를 이용하여 칼만필터 방법으로 앙상블 시간눈금을 생성하는 방법에 대해 논한다. 특히 앙상블 시간의 주파수안정도가 앙상블에 포함된 개별 시계의 주파수안정도보다 월등하게 우수해진다는 것을 보여주고자 한다. 미래 Korea Positioning System (KPS)의 시스템 시간을 위해 UTC(KRIS)가 안정적인 기준시간을 제공할 수 있도록 앙상블 시간눈금을 생성하는 것이 본 연구가 추구하는 최종목적이다. 단, 본 논문에서는 칼만필터 알고리즘과 그것을 구현하는 프로그램을 검증하는 차원에서 앙상블에 필요한 최소 시계 수인 3대의 시계를 이용하여 앙상블 시간눈금을 생성하고 그 특성을 살펴본다.

본 논문의 2장에서는 계산에 사용된 수학적 도구에 대해 설명한다. 시계 모델에 대해 이해하고, 앙상블 시간눈금 생성 방정식에 대해 알아본다. 그리고 일반적인 칼만 필터 알고리즘을 실제로 3-시계 시스템에 응용하는 방법을 보인다. 3장에서는 수소메이저 3대 사이에서 측정한 시간차이 데이터를 이용하여 칼만필터로부터 각 시계의 위상 및 주파수를 구하는 과정과 결과에 대해 설명한다. 4장에서는 칼만필터의 최종 결과물인 가중평균된 시간, 즉 앙상블 시간에 대해 주파수안정도 측면에서 논의한 후 5장에서 결론을 맺는다.

2. MATHEMATICAL TOOLS

2.1 Clock Model

원자시계와 같은 고정밀 주파수발생기의 주파수 또는 위상의 신호와 잡음을 분석할 때 진폭잡음은 무시할 만큼 작다고 가정한다. 그래서 위상과 주파수를 주요 분석 대상으로 삼는다. 어떤 순간 t에서의 주파수는 식 (1)과 같이 무차원의 상대주파수로 표현할 수 있다.

\(y(t)=\frac{v(t)-v_{0}}{v_{0}}=\frac{1}{2 \pi v_{0}} \frac{d}{d t} \phi(t)\)       (1)

여기서 ν₀는 명목주파수, ν(t)는 시점 t에서의 실제 주파수(단위: 헤르츠), ϕ(t)는 위상(단위: 라디안)을 나타낸다. 위 식 우변항에서 ϕ(t)/2πν₀≡x(t)로 정의하고, x(t)를 초(기호: s) 단위를 갖는 ‘위상’이라고 부른다. 따라서 식 (1)은 식 (2)와 같이 쓸 수 있다.

\(y(t)=\frac{d}{d t} x(t)\)       (2)

어떤 시계가 미래에 어떤 시간을 나타낼 것인지 예측하기 위해서는 해당 시계의 상태를 나타내는 x(t)와 y(t)에 관한 정보가 필요하다. 단, 여기서 말하는 시간이란 원자시계와 같이 성능이 우수한 시계가 나타내는 시간을 말하는 것으로 모두 1초 이내의 시간, 구체적으로 말하면 하루에 1 μs를 벗어나지 않는 시간을 대상으로 한다.

시점 (t-τ)에서 시점 t로 시간 τ가 흐른 후 시계의 상태를 나타내는 식 (3)을 흔히 ‘시계 모델’이라고 부른다 (Tavella 2008, Levine 2012).

\(\begin{array}{l} x(t)=x(t-\tau)+y(t-\tau) \cdot \tau+0.5 \cdot d(t-\tau) \cdot \tau^{2}+\xi \\ y(t)=y(t-\tau)+d(t-\tau) \cdot \tau+\eta \\ d(t)=d(t-\tau)+\zeta \end{array}\)       (3)

여기서 x(t), y(t), d(t)는 각각 시점 t에서 시계의 위상, 주파수, 주 파수 표류(frequency drift)를 나타낸다. 그리고 ξ, η, ζ는 랜덤 (random) 잡음을 나타내는데, 각각의 평균은 0이고 일정 분산을 가지는 정규분포를 하며, 서로 상관관계가 없다고 가정한다. 이 잡음들의 단위는 각 항의 좌변과 동일한 [s], [s/s], [s/s²]이다.

시계가 위치한 환경 조건(예: 온도, 전압 등)이 일정하게 유지된다면 식 (3)의 위상, 주파수 및 주파수 표류는 거의 일정하게 변하기 때문에 기준시계와 비교 측정하여 어느 정도 정확히 알아낼 수 있다. 그래서 이것들을 결정형(deterministic) 성분이라고 부른다 (단, 주파수 표류의 경우에는 그 값을 알아내기까지 긴 시간이 필요하기 때문에 아주 안정적인 기준시계 또는 주파수원이 있어야 한다). 이에 비해 랜덤 잡음은 확률적(stochastic) 성분으로 그 값을 정확히 예측할 수 없다. 결국 원자시계와 같은 고정밀 주파수 발생기의 단기 미래 상태를 예측할 때 예측 정확도는 랜덤 잡음 성분에 의해 결정된다 (장기 미래 예측에는 주파수 표류에 대한 정확한 정보가 필요하다. 하지만 본 논문에서와 같이 최대 평균시간이 106초 (약 12일) 수준인 경우에는 주파수 표류는 일정하다고 가정한다).

잡음을 분석할 때 주파수 영역에서 스펙트럼 밀도(spectral density)로 분석하거나 시간 영역에서 주파수안정도로 분석한다. 본 논문에서는 시간 영역에서의 분석법, 즉 평균시간(또는 적분시간)에 따른 알란편차(Allan deviation: ADEV)의 기울기로써 분석하는 방법을 사용했다 (Allan 1987).

연속적으로 측정한 위상(시간차이) 데이터로부터 알란편차를 계산하는 공식은 잘 알려져 있다. 최근에는 알란편차의 신뢰성을 높이기 위해 측정 데이터를 겹쳐서 계산하는 방법, 즉 overlapping ADEV를 많이 사용한다 (Riely 2008). 본 논문에서 말하는 알란편차는 모두 이 방법으로 계산한 것이다. 식 (4)는 알란분산(AVAR)을 계산하는 공식으로 이것의 제곱근이 알란편차(ADEV)이다.

\(\sigma_{y}^{2}(\tau)=\frac{1}{2(N-2 m) \tau^{2}} \sum_{i=1}^{N-2 m}\left(x_{i+2 m}-2 x_{i+m}+x_{i}\right)^{2}\)       (4)

여기서 괄호 속의 x_i 등은 연속적으로 측정한 시간차이 데이터를 나타내고, N은 가장 짧은 시간간격(τ₀)으로 측정한 데이터의 총 수이다. m은 평균시간(τ)을 결정하는 정수로서, τ = mτ₀의 관계를 가진다. 본 논문에서는 GPS에서의 측정시간 간격과 동일하게 τ₀ = 900 s를 선택했다 (Suess et al. 2010). 

원자시계의 종류에 따라 해당 시계에서 우세한 잡음의 종류는 다르다. 예를 들면, 세슘원자시계의 경우에는 White Frequency Modulation (WFM) 잡음과 Flicker FM (FFM) 잡음이 우세하다. 이에 비해 수소메이저의 경우에는 WFM과 Random Walk FM(RWFM) 잡음이 우세하다. 그래서 수소메이저에 대한 이론적 모델을 만들 때 WFM과 RWFM 잡음만 고려한다 (Tavella 2008). 이런 잡음 성분에 의한 주파수안정도 그래프, 즉 평균시간(τ)에 대한 알란편차 (σy (τ)) 그래프를 log-log 스케일로 그릴 때 그 기울기는 WFM은 (-1/2), FFM은 0, 그리고 RWFM은 (+1/2)이다. 수소메이저의 일반적인 특성은 단기 주파수안정도는 아주 우수하지만 장기 안정도는 RWFM의 영향으로 좋지 않다. 수소메이저에서 WFM과 RWFM이 구분되는 평균시간은 개별 시계에 따라 다르지만 대략 하루(86,400초)이다. 이에 비해 세슘원자시계에서 WFM과 FFM이 나누어지는 시간은 시계의 성능과 환경 조건에 따라 다르지만, 대략 수일∼수십일이다.

수소메이저의 주파수안정도(알란편차)를 구하기 위해서는 그 보다 우수한 기준 주파수원이 필요하다. 하지만 그런 주파수원이 없기 때문에 성능이 서로 비슷한 3대의 수소메이저끼리 서로 비교 측정하여 그들 사이의 상대적인 알란편차를 구한 후, 그것으로부터 수소메이저 각각의 절대 알란편차를 계산해내는 방법을 사용한다. 이 방법을 “3-cornered hat” 이라고 부른다 (Gray & Allan 1974, Riley 2008). 3대의 수소메이저를 각각 a, b, c라고 하고, σ2ab를 b를 기준으로 구한 a의 알란분산이라고 하면 3대의 수소메이저 사이의 상대 알란분산은 식 (5)와 같이 절대 알란분산으로 표현할 수 있다 (여기서 절대 알란분산이란 기준 주파수원의 성능이 월등하게 우수하여 그것의 잡음이 측정 대상에 영향을 미치지 않을 때 구한 알란분산을 의미함). 단, 수소메이저들이 서로 독립적이어서 각 신호들 사이에 상관관계가 없다고 가정할 때 이 관계식을 사용할 수 있다.

\(\sigma_{a b}^{2}=\sigma_{a}^{2}+\sigma_{b}^{2} ; \sigma_{a c}^{2}=\sigma_{a}^{2}+\sigma_{c}^{2} ; \sigma_{b c}^{2}=\sigma_{b}^{2}+\sigma_{c}^{2}\)       (5)

이 식들로부터 구한 개별 수소메이저의 절대 알란분산은 식(6)과 같다.

\(\begin{aligned} \sigma_{a}^{2} &=\frac{1}{2}\left(\sigma_{a b}^{2}+\sigma_{a c}^{2}-\sigma_{b c}^{2}\right) \\ \sigma_{b}^{2} &=\frac{1}{2}\left(\sigma_{a b}^{2}+\sigma_{b c}^{2}-\sigma_{a c}^{2}\right) \\ \sigma_{c}^{2} &=\frac{1}{2}\left(\sigma_{a c}^{2}+\sigma_{b c}^{2}-\sigma_{a b}^{2}\right) \end{aligned}\)       (6)

알란분산은 (-) 값이 될 수 없지만 식 (6)에서 보는 것처럼 경우에 따라서는 괄호 속이 (-)가 될 수 있다. 이 마이너스 분산은 성능이 너무 다른 시계들에 대해서 3-cornered hat을 적용하거나, 시계들 사이에 강한 상관관계를 가지고 있을 때 주로 장기(long-term) 안정도에서 발생한다. 그리고 측정 데이터 수가 충분히 많지 않을 때도 생길 수 있다. 본 논문에서는 위 식을 이용하여 절대 알란분산을 구할 때 IEEE에서 보급하는 Stable32라는 프로그램을 사용했다 (IEEE Stable32 2020).

2.2 Clock Ensemble Algorithm

같은 종류 동일 모델의 원자시계라 하더라도 시계들의 표시값(clock reading) 사이에는 미세한 차이가 난다. 시계 사이의 시간 차이 측정에 널리 사용되고 있는 시간간격 계수기(time interval counter)나 Dual Mixer Time Difference (DMTD) 등은 그 시간분해능이 보통 수 피코초(1 ps=10^-12 s)에서 수십 피코초이다. 오늘 날에는 이 보다 더 정밀한 시간차이를 알아내는 것이 필요해지고 있다.

여러 대의 시계 사이의 시간차이를 측정하고, 그 값들을 적절한 수학적 방법으로 평균하면 개별 시계들의 시간보다 더 안정적인 시간을 생성할 수 있다. 이런 개념으로 만든 시간이 앙상블 시간이다. N개의 시계들로 앙상블 시계를 구성하고, 그 중 어떤 시계 i의 표시값을 h_i 라고 하자. 그런데 시계 표시값 hi 는 측 정할 수 있는 양이 아니다. 시계 사이의 시간차이가 측정가능한(observable) 양이다. 이를 반영하여 앙상블 시간 TA를 어떤 시계 j의 표시값과의 차이로 나타내면 식 (7)과 같이 쓸 수 있다(Tavella 2008).

\(T A-h_{j}(t)=\frac{\sum_{i=1}^{N} p_{i}\left[h_{i}(t)-h_{j}(t)\right]}{\sum_{i=1}^{N} p_{i}}\)       (7)

여기서 p_i는 시계 i의 가중치(weight)이고, TA는 가중평균시간(weighted mean time)을 의미한다. hj 는 앙상블을 구성하는 시계들 중에서 기준(reference)으로 선택한 시계의 표시값을 나타낸다. 개별 시계들의 가중치 p_i 는 식 (8)과 같이 해당 시계의 알란분산(σ_y² (τ))에 반비례하는 값을 선택한다. 그렇게 함으로써 앙상블시간의 안정도를 높일 수 있다.

\(p_{i} \sim \frac{1}{\sigma_{y i}^{2}(\tau)}\)       (8)

여기서 τ는 식 (4)의 평균시간이고, 또한 식 (3)의 경과시간과 동일하다.

식 (8)의 가중치를 모든 시계에 대해 더하면, 앙상블 시계의 알란분산(σ²_ye(τ))은 개별 시계의 알란분산보다 작아진다. 왜냐하면, 아무리 안정도가 나쁜 (알란분산이 큰) 시계라도 그것들이 더해지면 식 (9)에서 보는 것처럼 전체 알란분산은 줄어들기 때문이다.

\(\frac{1}{\sigma_{y e}^{2}(\tau)}=\sum_{i=1}^{N} \frac{1}{\sigma_{y i}^{2}(\tau)}\)       (9)

알상블 시계를 구성하는 어떤 시계가 고장나서 앙상블에서 제외시키거나 새로운 시계가 앙상블에 포함되는 경우에도 앙상블 시간이 연속성을 유지하려면 특별한 방법이 필요하다. 앙상블을 구성하는 시계의 개수가 바뀔 때 발생하는 앙상블 시계의 시간 또는 주파수의 불연속(step)을 피하기 위해 시계들의 거동을 미리 예측하여 식 (7)에 포함시켜서 계산한다. 식 (3)의 시계 모델에서 결정형 파라미터 값들이 적절히 구해졌다면 미래 시점에서 시계가 어떤 시간을 나타낼 것인지 어느 정도 예측할 수 있고, 그 값을 식 (7)에 포함시키는 것이다. 새 시계가 도입되는 경우에도 그 시계의 파라미터 값을 충분히 정확히 파악하고, 예측한 후에 앙상블에 포함시킨다. 시점 t에서 시계 i의 예측 표시값을 h'_i (t)라고 한다면 식 (10)과 같이 계산한다.

\(T A-h_{j}=\frac{\sum_{i=1}^{N} p_{i}\left[h_{i}(t)+h_{i}^{\prime}(t)-h_{j}(t)\right]}{\sum_{i=1}^{N} p_{i}}\)       (10)

식 (10)이 앙상블 시간 TA의 실질적인 정의이다. 결국, 앙상블 시간은 기준시계와의 시간차이로 정의된다. 이 식에서 보는 것처럼 앙상블을 구성하는 개별 시계들의 가중치 pi 와 예측값 h'_i(t)을 결정하는 것이 앙상블 시간 생성에서 가장 중요하다. 가중치는 앞에서 설명한 것처럼 각 시계의 알란분산의 역수에 비례하도록 정한다. 한편 칼만필터 알고리즘은 바로 이 예측값을 결정하기 위해서 사용된다. 앙상블 시간을 생성한 후에 이것의 시간이나 주파수를 UTC(lab)이나 1차 주파수표준기에 맞추는(steering) 것은 또다른 일로서, 앙상블 시간 생성 알고리즘 밖에서 수행한다. 이에 대해서는 4장에서 다시 설명한다.

개별 가중치 p_i를 가중치들의 합으로 나누어 정규화해서 그 합이 1이 되도록 만들어 식 (10)을 다시 쓰면 식 (11)과 같다.

\(T A-h_{j}=\sum_{i=1}^{N} w_{i}\left[h_{i}(t)+h_{i}^{\prime}(t)-h_{j}(t)\right]\)       (11)

여기서 w_i는 각 시계의 상대 가중치를 나타내며, \(\sum_{i=1}^{N} w_{i}=1\)이다.

식 (11)에서 TA-h_j ≡x_j로 정의하고, 두 시계 표시값의 차이를 x_ij≡h_j -h_i 라 표시하며, 예측값 h'_i를 뒤에서 나올 칼만필터의 출력 (최종 추정값)을 나타내는 기호 \(\hat{x}_{i}\)로 표시하면 식 (11)은 식 (12)와 같이 쓸 수 있다.

\(x_{j}(t) \equiv T A-h_{j}=\sum_{i=1}^{N} w_{i}\left[\hat{x}_{i}(t)-x_{i j}(t)\right]\)       (12)

식 (10)이나 식 (12)는 앙상블 시간을 구할 때 반드시 필요하기 때문에 Basic Time Scale Equation (BTSE)이라고 부른다 (Tavella 2008).

2.3 Application of a Kalman Filter Algorithm to a 3-Clock System

칼만필터를 앙상블 시간눈금에 응용하기 위해서는 앞에서 살펴본 시계 모델 및 시간눈금 생성 알고리즘을 칼만필터의 해당 요소에 적절히 연결시키는 것이 필요하다. 칼만필터를 동작하기 위해 기본적으로 필요한 정보는 다음 4가지이다.

① 물리적 시스템의 전개 모델: 본 논문의 시계 모델에 해당

② 시스템 모델에 추가되는 잡음에 관한 정보: 잡음의 공분산 행렬로 표현됨

③ 초기값을 포함한 이전 추정값에 대한 정보: 상태 및 공분산 행렬의 추정값

④ 물리적 시스템과 연관된 측정값: 본 논문에서는 시계 사이 의 시간차이 측정값

물리적 시스템(예: 3-시계 시스템)의 어떤 상태(예: 시계의 위상)가 일정 시간 후에 어떻게 변했을 지 칼만필터로 추정해내려면 다음 과정을 반복한다: ①의 전개 모델에 의해 예측하고 추정했던 값과 ④의 현재 측정값이 반영된 값(구체적으로, Innovation)에 칼만이득을 적용하여 최종 추정값을 정한다. 이 최종 추정값이 다시 이전 추정값이 되어 같은 과정을 반복하게 되는데, 측정값에는 측정잡음이 추가되고, 추정값에는 그것의 불확실도(uncertainty)에 해당하는 공분산 행렬이 시간에 따라 진화하면서 따라다닌다. 이 반복 과정에서 이전 추정값과 현재 측정값이 가장 중요한 요소이다. 하지만 칼만필터의 성능은 시스템 잡음을 나타내는 공분산 행렬의 요소가 실제 시스템의 상황을 얼마나 정확히 나타내고 있는지에 의해 결정된다. 이런 요소들은 칼만필터를 동작시키기 전에 미리 결정해야 한다.

본 논문에서 사용하는 칼만필터 변수들의 기호와 의미가 Table 1에 정리되어 있다. 벡터와 행렬의 크기를 나타내는 n(=3N)과 m(=N-1)은, 3-시계 시스템인 경우에는 n = 9, m = 2이다. 단, N은 시계의 수이며, 본 논문에서는 N=3이다.

Table 1. The list of variables used in the Kalman filter of this paper. n=3N and m=N-1, where N is the number of clocks.

Symbol	Meaning	Size
xk	the state of the system at tk	(n×1) vector
\(\hat{x}_{k}\)	the estimate of xk	(n×1) vector
\(\hat{x}_{k}^{-}\)	the prediction of xk	(n×1) vector
zk	the measurement values at tk	(m×1) vector
wk	the system noise at tk	(n×1) vector
vk	the measurement noise at tk	(m×1) vector
Φ	the transition matrix	(n×n) matrix
H	the measurement matrix	(m×n) matrix
Γk	the covariance matrix of \(\hat{x}_{k}\)	(n×n) matrix
\(\Gamma_{k}^{-}\)	the prediction of Γk	(n×n) matrix
Q	the covariance matrix of wk	(n×n) matrix
R	the covariance matrix of vk	(m×m) matrix

 

표에 나타난 기호를 이용하여 칼만필터의 반복 과정을 설명하면 다음과 같다. 칼만필터는 tk-1과 tk로 표시한 두 이산(discrete)시점과 그 사이 시간 τ (단, τ = tk-tk-1)에서 반복적으로 동작한다. 즉 상태 예측(x̂ - k )과 이 상태와 관련된 실제 측정(zk), 그리고 이 두 가지 정보를 이용한 최종 상태추정값(x̂ _k)이다. 여기서 아래첨자 k 는 시점 t_k를 나타내고, 위 첨자(-)는 예측값을, 그리고 햇(^)은 추정값을 나타낸다. H는 측정 행렬을 의미하며, 예측값에 작용하여 이론적인 측정값을 구하는 데 사용된다. 결국 모델에 의해 이론적으로 예측한 측정값과 실제 측정값의 차이가 최종 추정값을 결정하는 데 기여한다.

칼만필터의 계산 순서는 Fig. 1에 나와 있다. 여기서 시작점 t₀에서의 상태 추정값 x̂ ₀ 및 공분산 추정값 Γ₀는 다른 방법으로 미리 구해서 주어져야 한다. 그림에 나와 있는 공식과 순서대로 계산이 진행되는데, 3번째 단계에서 측정이 수행되어 입력되고, 해당 사이클의 최종 추정값이 출력으로 나온다.

Fig. 1. Four steps and equations of the recursive Kalman filter algorithm.

앙상블 시계를 구성할 때 최소 시계의 개수는 N=3이다. 본 논문에서는 세 대의 수소메이저를 사용했는데, 각각 1, 2, 3번으로 구분하고 1번 시계를 기준시계로 사용한다. 개별 시계의 상태를 식 (3)의 x, y, d로 나타낸다면 이들로 구성된 시스템의 상태 벡터S는 식 (13)과 같이 9×1의 열벡터(column vector)로 표현할 수 있다 (상태 벡터를 지금까지 x_k로 표시했지만 시계의 위상을 나타내는 x와 혼동을 피하기 위해 여기서는 S로 표시함).

\(S=\left[x_{1}, y_{1}, d_{1}, x_{2}, y_{2}, d_{2}, x_{3}, y_{3}, d_{3}\right]^{T}\)       (13)

시계 번호 (1, 2, 3)에 따라 시계 상태(x, y, d)를 순서대로 나열했다. 이 순서에 맞추어 뒤에 나오는 여러 행렬에서 요소들의 순서도 정해진다.

시스템의 전이 행렬 Φ는 시점 t_k-1에서 시점 t_k로 시스템의 상태 변화(진화)를 나타내는 것으로 결국 식 (3)의 시계 모델을 행렬로 표현한 것이다. 시계 하나 당 3개의 방정식이 필요하므로 Φ는 9×9의 정방 행렬이다. 이 행렬은 개별 시계 i에 대한 3×3 행렬Φ_i (i=1, 2, 3)가 식 (14)와 같이 대각 행렬을 이룬다. 단, 오른쪽 식 에서 0 하나는 0으로 구성된 (3x3)의 정방 행렬이다.

\(\boldsymbol{\Phi}_{i}=\left[\begin{array}{ccc} 1 & \tau & \tau^{2} / 2 \\ 0 & 1 & \tau \\ 0 & 0 & 1 \end{array}\right], \quad \boldsymbol{\Phi}=\left[\begin{array}{ccc} \boldsymbol{\Phi}_{1} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \boldsymbol{\Phi}_{2} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \boldsymbol{\Phi}_{3} \end{array}\right]\)       (14)

측정 행렬 H는 상태 벡터 S에 작용하여, 측정한 시계 사이의 시간차이를 나타낸다. 따라서 식 (15)와 같이 2×9의 행렬로 표현된다.

\(\boldsymbol{H}=\left[\begin{array}{ccccccccc} 1 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 \end{array}\right]\)       (15)

시계 1번을 기준시계로 선택했으므로, 시계들 사이의 시간차 이 측정값은 x₁-x₂, x₁-x₃이다. 따라서 HS는 식 (16)과 같이 2×1의 벡터이다.

\(\boldsymbol{H} \boldsymbol{S}=\left[\begin{array}{l} x_{1}-x_{2} \\ x_{1}-x_{3} \end{array}\right] \equiv\left[\begin{array}{l} x_{21} \\ x_{31} \end{array}\right]\)       (16)

여기서 우변 마지막 항은 기호를 단순하게 표현하기 위해 정의한 것으로, 식 (12)의 x_ij ≡h_j-h_i와 연결돤다. 그리고  x₂₁=-x₁₂, x₃₁=-x₁₃의 관계가 성립한다.

시스템 잡음의 공분산 행렬은 \(\boldsymbol{Q}=E\left[\boldsymbol{w}_{k} \boldsymbol{w}_{k}^{T}\right]\)로 표현되고, 잡음 w_k는 식 (3)의 잡음 성분 ξ, η, ζ의 조합으로 구성된다. 그런데 이 잡음들은 시간이 경과함에 따라 확산되는데, 잡음의 종류에 따라 확산계수(diffusion coefficient)가 다르다. 확산계수로써 공분산 행렬 Q를 표현한 식은 Brown (1991)과 Hutsell (1995)에 나와 있다. 시계 i의 공분산 행렬은 3×3의 Qi (단, i=1, 2, 3)로 표현하고, Qi 의 구성요소는 확산계수 q_wf , q_rw, q_rr로 표현한다. 여기서 q_wf는 White FM, q_rw는 Random Walk FM, q_rr는 Random Run FM 잡음의 확산계수이다. 시스템의 공분산 행렬 Q는 식 (17)과 같은 Q_i 3개로 구성된, 9×9의 정방 행렬이다.

\(\boldsymbol{Q}_{i}=\left[\begin{array}{ccc} q_{w f}+q_{r w} \tau^{3} / 3+q_{r r} \tau^{5} / 20 & q_{r w} \tau^{2} / 2+q_{r r} \tau^{4} / 8 & q_{r r} \tau^{3} / 6 \\ q_{r w} \tau^{2} / 2+q_{r r} \tau^{4} / 8 & q_{r w} \tau+q_{r r} \tau^{3} / 3 & q_{r r} \tau^{2} / 2 \\ q_{r r} \tau^{3} / 6 & q_{r r} \tau^{2} / 2 & q_{r r} \tau \end{array}\right]\)       (17)

같은 종류의 잡음이라도 시계들마다 q_wf , q_rw, q_rr는 서로 다른 값을 가진다. 따라서 본 연구에서는 총 9개의 q 값을 알아야 하는데, 이미 발표된 값을 Table 2에 나타냈다 (Suess et al. 2010). 표에서 Active H-maser는 본 논문에서 사용한 능동형 수소메이저와 같은 종류이다. 세슘원자시계와 비교해보면 q_wf는 3 차수 (order) 작아서 단기 안정도가 아주 우수하다는 것을 짐작할 수 있다. 이에 비해 장기 안정도와 관련된 q_rw와 q_rr은 오히려 세슘원자시계보다 2 차수 더 크다. 즉, 장기 안정도는 세슘원자시계가 더 우수하다. 분수시계(Fountain clock)의 확산계수가 다른 시계보다 훨씬 작기 때문에 장기와 단기 안정도 모두 월등히 우수하다. 이에 비해 GPS, Galileo, BeiDou 위성에 탑재되어 있는 루비듐원자시계는 단기 안정도는 세슘원자시계보다 좋으나 장기 안정도는 시계들 중에서 가장 나쁘다 (Galileo와 BeiDou에는 수동형 수소메이저도 같이 탑재됨).

Table 2. Typical values of components of the system noise covariance matrix (Suess et al. 2010).

Clock	qwf (unit: s)	qrw (unit: 1/s)	qrr (unit: 1/s3)
Cs atomic clock	2.5×10-23	4.4×10-37	5.0×10-53
Active H-maser	2.8×10-26	1.1×10-35	4.4×10-51
Fountain clock	2.5×10-26	1.1×10-37	1.1×10-55
Rb atomic clock	1.0×10-24	1.1×10-35	2.8×10-46

 

시계의 잡음별 확산계수가 시계의 주파수 안정도에 영향을 미치는데, 이 q 값들과 알란분산(AVAR)은 식 (18)과 같은 관계가 있다 (Chaffee 1987, Hutsell 1995).

\(\sigma_{y}^{2}(\tau)=\frac{q_{w f}}{\tau}+\frac{q_{r w} \tau}{3}+\frac{q_{r r} \tau^{3}}{20}\)       (18)

Greenhall (2003)은 식(18)를 이용하여 칼만필터 알고리즘 안에서 알란분산을 구하고, 그로부터 각 시계의 가중치를 결정했다. 그런데 이 q 값들은 시스템이 바뀌지 않는 한 동일한 값을 가지므로 알란분산과 가중치도 변하지 않는다. Galleani & Tavella (2003)는 가중치가 변하지 않는다면 식 (11)의 예측값 h'_i(t)는 아무 쓸모가 없다고 주장한다. 세 잡음의 확산계수 중 q_rr은 본 논문의 관찰 영역인 평균시간 10³초에서 10⁶초 사이의 알란편차에는 영향을 미치지 않는다. 결국 WFM과 RWFM 잡음에 의해서만 알란편차 그래프가 결정된다. 그렇기 때문에 식 (18)로 구한 알란분산으로 시계의 가중치를 결정하는 것은 시계의 실제 상황을 반영하지 못할 수 있다. 그래서 본 논문에서는 칼만필터로 구한 각 시계들의 위상(시간차이) 값과 식 (4)를 이용하여 알란분산을 계산했다.

한편, 식 (12)에서 N = 3일 때 식을 전개해보면 식 (19)와 같이된다. 단, 1번 시계가 기준시계이다.

\(\begin{aligned} x_{1}(t) &=T A(t)-h_{1}(t)=\sum_{i=1}^{3} w_{i}\left[\hat{x}_{i}(t)-x_{1 i}(t)\right] \\ &=w_{1}\left[\hat{x}_{1}-x_{11}\right]+w_{2}\left[\hat{x}_{2}-x_{21}\right]+w_{3}\left[\hat{x}_{3}-x_{31}\right] \\ &=w_{1} \hat{x}_{1}+w_{2} \hat{x}_{2}+w_{3} \hat{x}_{3}-w_{2} x_{21}-w_{3} x_{31} \end{aligned}\)       (19)

여기서 x₁₁ =0이고, x₂₁과 x₃₁은 식 (16)처럼 각각 시계 2와 1, 시계 3과 1 사이의 시간차이를 측정한 값이다 (이 식에서 x₂₁과 x₃₁로 표시한 측정값을 추정값 등과 혼동을 피하기 위해 Fig. 1의 칼만필터 알고리즘에서는 z_k로 표시했음). w₁ , w₂, w₃은 각각 시계 1, 2, 3의 알란분산(AVAR)의 역수에 비례하는 값으로 결정한다. 그리고 x̂₁, x̂₂, x̂₃은 칼만필터에서 각 사이클마다 최종 추정한 시계 1, 2, 3의 위상값이다. 결론적으로, 앙상블 시간 TA는 기준시계 1번의 표시값과의 차이인 x₁(t)에 의해 정의되고, 이것은 칼만필터로 추정한 세 시계의 위상값과 시계 사이의 시간차이 측정값에 대한 가중평균으로 결정된다.

3. DATA ACQUSITION AND PROCESSING

KRISS에서 UTC(KRIS)를 생성하는 과정과 수소메이저에 대해 소개한 바 있다 (Lee et al. 2020). 본 논문에서 사용한 3대의 수소메이저는 각각 H8, H6, H4로 명명했다. 이 중 2019년 8월 29일(MJD 58724)부터 KRISS의 마스터 시계로 사용하고 있는 H8을 본 논문에서 기준시계로 사용했다. UTC(KRIS)를 생성하고, 시계들을 비교 측정하기 위한 시스템 개략도는 Fig. 2와 같다. 마스터 시계에서 나오는 5 MHz 주파수를 기준으로 다른 시계들에서 나오는 5 MHz 주파수를 Multi-channel Measurement System (MMS)에서 1분 마다 비교 측정하여 그 시간차이 데이터가 PC에 저장된다. 그 중에서 필요한 데이터를 다운로드하여 계산에 사용했다.

Fig. 2. Schemetic diagram of the system for time comparison, data acquisition, and generation of UTC(KRIS).

시계 비교 데이터는 MJD 58779(2019.10.23)부터 58923(2020.3.15)까지 145일 간의 데이터를 사용하되 15분(900초)마다 측정한 것을 추출했다. 세 시계 사이의 시간차이, 즉 (H6-H8), (H4-H8), (H4-H6)의 데이터들이 시작점(MJD 58779.000)에서 0이 되도록 나노초(ns) 단위의 값에서 각각 정수 부분을 빼거나 더 했다. 이렇게 구한 시간차이 데이터 13,920개 × 3 세트를 그린 것이 Fig. 3이다. 145일 동안 대략 1 μs∼2.5 μs 만큼 변한 것을 알 수 있다. 이것을 이 기간동안의 상대주파수로 나타내면 약 0.8×10^-13 ∼ 2×10^-13 에 해당한다.

Fig. 3. Time differences of 3 pairs of hydrogen masers for 145 days from MJD 58779 (Oct. 23, 2019) to MJD 58923 (Mar. 15, 2020).

이 시간차이 데이터를 알란분산을 구하는 식 (4)에 대입하여 구한 알란편차(ADEV)가 Fig. 4이다. 그림은 최소 평균시간 900초부터 최대 약 10⁶초까지 나타내었다. (H4-H8)의 단기 안정도가 우수하고, 장기 안정도는 (H6-H8)이 우수함을 알 수 있다. H8이 포함된 그래프의 안정도가 (H4-H6)에 비해 더 좋게 나왔다. 즉, H8이 다른 두 시계보다 안정도가 우수하다는 것을 짐작할 수 있다.

Fig. 4. Allan deviation curves of 3 pairs of clocks, calculated from the measured time difference shown in Fig. 3.

세 시계 사이의 알란편차 그래프로부터 식 (6)을 이용하여 3-cornered hat 방법으로 시계 각각의 절대 주파수안정도를 구했다. 그런데 일부 시계는 평균시간 10⁶초 부근에서 마이너스 분산이 나왔다. 그래서 1분마다 측정한 145일간의 데이터, 즉 263,519개 × 3 세트의 데이터를 이용하여 알란편차를 다시 계산했고, 그 중 비교 구간의 것만 나타낸 것이 Fig. 5이다. 그림에서 보듯이 H8이 10^-16 수준의 단기 및 장기 안정도를 보이며 세 시계 중 가장 우수하다. 이 그래프들은 뒤에서 칼만필터의 출력인 세 시계의 위상값으로부터 구한 알란편차와 비교할 때 기준으로 사용된다.

Fig. 5. Allan deviation curves obtained by a 3-cornered hat method from the Allan deviations of Fig. 4.

본 논문에서 칼만필터와 알란편차 등은 모두 MatLab으로 프로그래밍하여 계산했다.

칼만필터 알고리즘(Fig. 1)에서 입력되는 측정값 z_k는 Fig. 3의 시간차이 데이터 중 (H6-H8)과 (H4-H8)이다. 그리고 이 값들은 식 (16)과 (19)에서 x₂₁과 x₃₁에 해당한다. 즉, H8이 1번 기준시계이고, H6이 2번 시계, H4가 3번 시계이다.

3-시계 시스템에 대한 칼만필터 계산을 수행하기 위해서 값이 미리 주어져야 하는 파라미터들은 28개인데, 그것들은 다음과 같이 4종류로 분류할 수 있다.

① 시스템 잡음 공분산 행렬 Q의 구성 요소: Table 2 참조

-  시계 3대의 잡음별 확산계수: 총 9개

② 측정 잡음 공분산 행렬 R의 구성 요소

- 측정 시스템을 고려한 잡음의 분산: 총 1개

③ 시스템의 상태 벡터 S의 초기값: 식 (13) 참조

- 시계 3대의 위상(x), 주파수(y), 주파수 표류(d)의 초기값: 총 9개

④ 상태 추정값의 공분산 행렬 Γ의 초기값

- 시계 3대의 위상(x), 주파수(y), 주파수 표류(d) 추정값의 분산값: 총 9개

이 중 ③번은 해당 시계의 지난 이력이나 현재 상태로부터 그 값을 어느 정도 추정할 수 있다. 시계 위상의 경우는 데이터를 준비하는 과정에서 시작점이 나노초 단위에서 정수 부분이 0이 되도록 만들었기 때문에 쉽게 정할 수 있다. 주파수 표류는 추정값에서 한 차수를 변화시켜도 평균시간 τ∼10⁶ s에서 미미한 변화만 나타났다. 이에 비해 주파수에 따른 변화는 크게 나타났다. 하지만 초기값의 불일치에 의한 영향은 칼만필터가 처음 동작할 때 크게 나타나고 반복 수행하면서 금방 사라진다. 자세한 것은 다음 장에서 설명한다.

④번은 ③번 추정값의 분산에 해당하는 것으로 그 초기값은 아주 작다. 대략 추정한 분산값에서 한 두 차수 변화시켜도 칼만필터의 출력(예: 시계의 위상 추정값)에는 별 영향이 없었다. 하지만 이 값은 칼만필터가 반복 수행할수록 증가하게 된다.

②번은 측정 시스템에 따라 전혀 다른 값을 가질 수 있다. Fig. 2에 표시된 MMS는 기본적으로 Dual-Mixer Time Difference 방법으로 시간차이를 정밀하게 측정한다. 이 방법으로 인해 측정에 참여한 한쌍의 시계 사이에서 상관관계가 발생한다는 것은 이미 발표한 바 있다 (Lee et al. 2020). 우리와 동일한 방법을 사용한 다른 기관의 결과를 본 논문에서 그대로 사용했지만 (Levine 2012), 그 값을 한 두 차수 바꾸더라도 결과에는 별 차이가 없었다.

값을 결정하는 데 제일 어려운 것은 ①번이었다. 본 논문에서는 알란편차 그래프 모양의 변화를 관찰하면서 이 확산계수들의 값을 결정했다. 구체적으로 말하면, WFM과 RWFM 및 RRFM의 확산계수인 q_rw, q_wf , q_rr의 값에 따라 칼만필터로부터 시계의 위상 추정값을 구하고, 이 값을 이용하여 알란편차 그래프를 그렸다. 이때 Fig. 5에 나와 있는 각 시계들의 절대 주파수안정도(ADEV)그래프를 기준으로 삼았다. 즉 알란편차 그래프가 Fig. 5와 비슷한 모양이 되었을 때의 값을 찾았다. 그런데 시계들의 q 값은 상호 영향을 미쳤다. 예를 들면, 어떤 시계의 q 값을 변화시키면 다른 시계의 알란편차 그래프에서 변화가 나타나는 등, 상대적인 변화를 보였다. 이 q 값들은 시스템을 구성하는 시계가 바뀌거나 고장이 나서 고치거나 하지 않는 이상 한번 그 값이 결정되면 계속 사용한다.

칼만필터가 동작하는 데 필요한 파라미터의 값들이 시스템의 상태를 얼마나 정확히 나타내느냐에 따라 칼만필터의 유용성과 효과성이 결정된다. 그래서 본 논문에서는 일차적으로 이 파라미터 값을 가능하면 정확히 결정하는 데 주의를 기울였다. 그리고 칼만필터는 추후 실시간(real-time)으로 앙상블 시간을 생성하는데 사용될 것인데 그러기 위해서는 일정 주기로 알란분산을 계산하고 그것에 반비례하는 가중치를 일정 주기로 업데이트 해야 한다. 하지만 본 논문에서는 우선 145일간 전체 데이터를 이용하여 알란분산을 구하고, 그것으로부터 가중치를 계산했다. 그리고 그 가중치를 미세 조정하면서 앙상블 시간의 주파수안정도를 비교 평가했다.

4. RESULTS AND DISCUSSION

앞에서 설명한 28개의 칼만필터 파라미터에 대한 최적 값을 결정했다. 그 중 시스템 잡음 공분산 행렬의 구성 요소 q 값과 시스템의 상태 x, y, d의 초기값을 정리한 것이 Table 3이다. 이 값들을 바탕으로 칼만필터로써 각 시계들의 위상값을 구한 후 그것으로부터 알란편차를 계산했다. q 값 중에서 WFM과 RWFM의 확산계수(q_wf, q_rw)에 대해 알란편차 그래프는 민감하게 변했다. 이에 반해 RRFM의 확산계수(q_rr)는 Table 3의 값에서 한 두 차수 변화시켜도 그래프에서는 평균시간 τ∼10⁶ s 부근에서 미미한 변화만 있었다. 하지만 이보다 장기 안정도에서는 그 영향이 일정 크기로 나타날 것으로 예상한다. Fig. 5와 가장 비슷한 모양을 구한 결과가 Fig. 6이다.

Table 3. Values of the system noise covariance matrix components (qs) and initial values of the clock states (x, y, and d) used in this Kalman filter.

Clock No.	qwf	qrw	qrr	x	y	d
Clock 1 (H8)	1.0×10-26	5.5×10-35	3.0×10-51	1.0×10-10	-1.5×10-13	0.5×10-22
Clock 2 (H6)	5.0×10-26	8.0×10-35	5.0×10-51	1.0×10-10	-4.5×10-14	1.0×10-21
Clock 3 (H4)	3.0×10-26	9.5×10-35	2.0×10-51	1.0×10-10	9.0×10-14	1.0×10-22

 

Fig. 6. Allan deviation curves of the 3 clocks, obtained from the phases estimated by the Kalman filter.

칼만필터 출력물로서 위상 대신에 주파수 y를 선택하여 그린 그림이 Fig. 7이다. Y-축의 눈금은 무차원의 10^-13으로 위상에 비해 훨씬 정밀하게 변화를 관찰할 수 있다. 이 그림은 변동이 큰 초기 데이터 29개를 삭제한 것이지만 MJD 0 부근에서는 아직 그 잔재가 있는 것을 볼 수 있다. H6은 처음에 아래로 값이 툭 떨어졌다가 안정 상태가 되었고, H8은 위로 갑자기 솟았다가 안정 상태가 된 것을 볼 수 있다. MJD 135일 부근에서 3개 그래프 모두 작은 피크가 있는 것을 볼 수 있다. 세 시계에서 동시에 튀었다는 것은 기준 시계(H8)에서 어떤 변화가 있었다는 것을 의미한다.

Fig. 7. The frequency y of each clock estimated by the Kalman filter.

칼만필터 반복 과정의 매 사이클에서 시스템 상태의 최종 추정값을 결정하는 데 사용되는 Innovation은 측정값의 변화를 더욱 정밀하게 보여준다. Fig. 8은 두 측정값 (z86 및 z84)에 대한 Innovation 그래프이다. 두 그래프를 구분해서 볼 수 있도록 z86그래프를 +2 × 10^-11 만큼 위로 이동시켰다. 그림에서 보는 것처럼 Innovation은 거의 백색잡음의 형태를 나타낸다. MJD 0 부근에서 z86 그래프는 급격히 감소하는 데, Fig. 7과 마찬가지로 초기값의 불일치에 의한 영향이 남아있기 때문이다. 그런데 MJD 135 부근에서 두 그래프 모두 크게 튄 것을 볼 수 있다. 주파수 그래프(Fig. 7)에서 보이는 작은 피크가 여기서는 훨씬 크게 보이는 것이다. 또한 142일 부근에서는 Fig. 7에서는 보이지 않던 튐을 볼 수 있다. 이처럼 Innovation을 관찰하면 측정값에서 발생하는 이상값(outlier)를 쉽게 찾을 수 있다. 이것은 향후 실시간 시간눈금을 생성할 때 이상값 검출에 활용할 수 있을 것이다.

Fig. 8. Innovations with two measured time difference data, z86 and z84. The Innov-z86 is shifted by +2×10^-11 for clarity

3-시계 시스템에서 1번 시계를 기준으로 앙상블 시간을 구하는 식이 식 (19)이다. 이 식은 각 시계의 위상 추정값과 시간차이 측정값에 가중치를 부여함으로써 앙상블 시간을 구하는 것을 나타낸다. 가중치는 식 (8)과 (9)에서 설명한 것처럼 각 시계의 알란분산에 반비례하는 값으로 결정한다. 본 논문에서는 Fig. 6에 나타난 알란편차 그래프에서 평균시간 57,600초에서 알란분산(AVAR) 분의 1의 값, 즉 식 (8)의 p_i를 각각 구하고 세 시계의 가중치의 합계(p₁+p₂+p₃)를 구했다. 이 합계로써 각각의 가중치 p_i를 나누어 상대 가중치 wi 를 결정했다. 여기서 57,600초는 세 시계의 알란편차가 대략 최소가 되는 평균시간을 선택한 것이다. 이렇게 구한 가중치는 w1=0.6055, w2=0.1876, w3=0.2069 이었다. 이 가중치를 식 (19)에 대입하고, 1번 시계를 기준으로 가중평균된 시간 차이 x1 (t)를 구한 결과가 Fig. 9의 WMH8로 표시된 그래프이다. 나머지 3개의 그래프는 비교를 위해 같이 그렸다.

Fig. 9. Phase of the three member clocks and their weighted-mean ensemble clock with the reference clock of H8. For the WMH8 curve, the weights are w_H8=0.6055, w_H6=0.1876, and w_H4=0.2069, respectively

가중평균된 시간차이 데이터 WMH8를 식 (4)에 대입하여 계산한 알란편차 그래프가 Fig. 10의 맨 아래 그래프이다. 나머지 3개 그래프는 Fig. 6과 동일하며 비교를 위해 같이 그렸다. 그림에서 보는 것처럼 WMH8의 알란편차가 세 시계 각각의 알란편차보다 작은데, 평균시간 57,600초에서 최소값이 1.2×10^-16에 이른다. 앞에서 구한 가중치를 보면 1번 시계 (H8)의 가중치가 전체의 60%를 넘는다. 나머지 두 시계(H6, H4)는 각각 19%와 21%를 차지하여 거의 비슷한 비율이다. 주파수안정도가 가장 좋은 H8에 대한 의존성이 가장 높게 나온 것은 당연하다. 하지만 이 비율이 앙상블 시간의 안정도 측면에서 최적인지 확인하기 위하여 가중치 값을 미세 조정하면서 앙상블의 알란편차가 최소가 되는 값을 찾아보았다.

Fig. 10. Allan deviation of the weighted-mean ensemble clock calculated with the phase of Fig. 9. The minimum ADEV of the WHM8 curve is 1.2×10^-16 at the averaging time of 57,600 s.

앙상블 시계의 알란편차는 Fig. 10에서 보는 것처럼 평균시간 (τ)이 불연속적이다. 따라서 평균시간 57,600초에서의 알란편차가 최소값이 아닐 가능성이 있기 때문에 가중치 값을 조금씩 바꾸어 가면서 알란편차 변화량을 관찰하였다. 이 과정에서 특정 가중치에서 거의 한차수 작은 알란편차를 나타낸다는 것을 발견했다. 즉 가중치의 값이 각각 w₁=0.7206, w₂=0.1452, w₃=0.1342 일 때 Fig. 11과 같이 최소 알란편차는 약 2×10^-17로서 Fig. 10의 결과에 비해 6분의1에 해당하는 작은 값을 나타냈다. 이것의 시간 안정도 (TDEV) σ_x(τ)=τσ_y(τ)/√_3 는 하루 동안에 수 피코(10^-12) 초 이내에서 변한다는 것을 의미한다.

Fig. 11. Allan deviation of the weighted-mean ensemble clock when the weights of three member clocks are w_H8=0.7206, w_H6=0.1452, and w_H4=0.1342, respectively. The minimum ADEV is 2×10^-17.

Fig. 11의 결과는 처음 예상했던 것보다 훨씬 작은 값으로, 이렇게 되는 이유를 아직 찾지 못했다. 실험적으로 이 결과가 맞는지를 확인하기 위해 현재 운영되고 있는 측정 및 조정 시스템(Fig. 2)과 비슷한 시스템을 구성하여 앞으로 실험할 예정이다. 이와 더불어 이론적으로 그 이유를 밝히는 연구가 필요하다. 특히, 시계 모델에서 기본 가정으로 채택한 잡음성분들 사이의 무상관성(uncorrelation)은 실제 상황에서는 사실이 아니라는 것이 밝혀졌다 (Lee et al. 2020). 또한 상관성을 갖는 잡음들은 식 (6)으로 절대 알란편차를 계산할 때 실제 보다 낮은 값을 나타낼 수 있다. 이런 점들을 고려하여 잡음 사이의 상관관계가 알란편차에 미치는 영향에 대해서 심도 있는 연구가 더 필요하다.

이렇듯 시계들의 가중평균을 구하거나 가중치를 미세 조정함으로써 단기 주파수안정도가 개선된다는 것은 소프트웨어적으로 앙상블 시간을 생성하는 것이 아주 우수한 시계를 개발한 것과 동일한 효과를 볼 수 있다는 것을 뜻한다. 새로운 시계를 개발하는 데 많은 시간과 노력, 그리고 비용이 소요되는 것을 감안한다면 앙상블 시간을 생성하는 것은 경제적인 측면에서도 아주 바람직한 일이다. 그런데 앙상블을 구성하는 시계 중 어떤 하나의 시계에 가중치가 집중되는 것은 앙상블을 구성하는 의미가 줄어드는 것이다. 앙상블을 구성하는 시계의 개수가 더 늘어나는 경우에는 시계의 최대 가중치 값에 의도적인 제한을 두는, 이른바 관리적 한계(administrative limt)를 설정하는 것이 필요하다 (Levine 2012).

본 논문에서는 기준시계로서 단기 안정도가 우수한 H8을 선택했다. 하지만 다른 시계(H6 또는 H4)를 기준시계로 선택하더라도 가중평균된 시간은 Fig. 10과 동일한 결과가 나온다는 것을 확인했다. 다시 말하면, 기준시계가 되기 위한 특별한 조건이 없다는 것이고, 이것이 앙상블 시간생성 알고리즘의 큰 장점이다. 단지, 기준시계를 자주 바꾸는 것은 바람직하지 않기 때문에 수명이 많이 남아 있는 시계를 선택하는 것이 바람직하다.

주파수안정도가 우수한 가중평균시간(즉 앙상블 시간)을 생성한 다음에 할 일은 이 시간이 UTC와 일치되도록 만드는 것이다. 이 과정은 칼만 필터 알고리즘 밖에서 수행되는 것으로 외부로부터 기준값이 제공되어야 한다. BIPM은 UTC와 UTC(lab)의 시간차이를 매달 발표한다. 세계 각국의 시간표준연구실(lab)에서 한달 동안 측정한 UTC(lab)-clock 데이터를 BIPM에 보내면 BIPM에서는 15일 후에 UTC-UTC(lab)를 Circular T를 통해 발표한다(Circular-T 2020). lab 입장에서는 첫 번째 데이터를 측정한 지 45일만에 UTC-UTC(lab) 값을 알게 된다. 이 값을 알면 해당 연구실에서는 UTC-clock 값을 계산할 수 있다. H8을 clock으로 선택하고, BIPM으로부터 UTC-H8=δ를 얻었다고 하자. 식 (19)에서 TA-H8=WMH8을 구했으므로 이 두 식을 빼면 UTCTA=(δ-WMH8)이 된다. 만약 우변값을 0으로 만든다면 TA=UTC 가 된다. 즉 앙상블 시간 TA는 UTC를 따라가게 된다. 여기서 (δ-WMH8)을 0으로 만드는 장치로 Fig. 2의 Microphase stepper가 사용된다. 이 장치는 입력 주파수의 위상 또는 주파수를 미세하게 조정한 후 출력으로 내보내는 장치이다. 따라서 입력단에 기준시계 H8의 주파수를 입력하고, (δ-WMH8)가 0이 되도록 이 장치를 조정하면 그 출력은 UTC에 맞게 조정된 앙상블 시간이 나온다. 이 앙상블 시간을 UTC(KRIS)으로 사용하면 UTC에 근접하고, 안정도가 높은 앙상블 시간눈금이 생성된 것이다. 단, (δ-WMH8)를 조정할 때 한번에 0으로 만들면 그 순간에 주파수 또는 위상이 크게 변하여 안정도가 나빠지므로 긴 시간(예: 30일)에 걸쳐 나누어서 조정하는 것이 필요하다.

5. CONCLUSIONS

본 논문에서는 미래 KPS의 시스템 시간의 기준이 될 UTC(KRIS)를 더욱 안정적이고 강건하게 생성하는 방법에 대해 연구했다. KRISS에서 운용 중인 세 대의 수소메이저를 이용하여 칼만필터 기반의 앙상블 시간눈금을 생성하고 그 주파수안정도를 높이는 방법을 찾았다. 이를 위해 우선 3-시계 시스템의 여러 파라미터 값을 결정했고, 다음으로 시계들의 알란분산에 반비례하는 가중평균으로 만든 앙상블 시계의 주파수안정도(알란편차)가 개별 시계들의 주파수안정도보다 개선된다는 것을 확인했다. 특히 시계의 가중치를 단순히 알란분산에 반비례하는 값으로 결정하기보다 그 값 주변에서 미세 조정하면 훨씬 더 안정도가 개선된다는 것을 발견했다. 구체적으로 말하면, 일반적인 방법으로 구한 가중치를 사용하여 구한 앙상블 시계의 알란편차는 평균시간(τ) 57,600초에서 σ_y(τ)=1.2×10^-16이었지만, 미세 조정하여 찾은 가중치에서는 σ_y(τ)=2×10^-17으로, 약 6분의 1로 줄어들었다. 이것의 시간 안정도(TDEV)는 하루 동안에 수 피코(10^-12) 초에 해당한다. 하지만 이렇게 우수하게 보이는 결과가 잡음 사이의 상관관계에 의해서 나타난 현상일 수 있기 때문에 실험 및 이론으로 검증하는 연구가 필요하다. 향후 실시간 시간눈금을 실제로 생성할 때 이런 특성을 고려하여 설계하면 더욱 안정적인 위성항법용 기준시 생성 시스템을 구현 할 수 있을 것으로 기대된다.