Comparison of AT1- and Kalman Filter-Based Ensemble Time Scale Algorithms

Abstract

We compared two typical ensemble time scale algorithms; AT1 and Kalman filter. Four commercial atomic clocks composed of two hydrogen masers and two cesium atomic clocks provided measurement data to the algorithms. The allocation of relative weights to the clocks is important to generate a stable ensemble time. A 30 day-average-weight model, which was obtained from the average Allan variance of each clock, was applied to the AT1 algorithm. For the reduced Kalman filter (Kred) algorithm, we gave the same weights to the two hydrogen masers. We also compared the frequency stabilities of the outcome from the algorithms when the frequency offsets and/or the frequency drift offsets estimated by the algorithms were corrected or not corrected by the KRISS-made primary frequency standard, KRISS-F1. We found that the Kred algorithm is more effective to generate a stable ensemble time scale in the long-term, and the algorithm also generates much enhanced short-term stability when the frequency offset is used for the calculation of the Allan deviation instead of the phase offset.

1. INTRODUCTION

위성항법시스템에서 가장 기본적이고 중요한 요소는 정확한 시간을 생성하여 방송하는 것이다. 이를 위해 항법 위성에는 원자시계들이 탑재되고, 이 시계들을 모니터링하고 조정하는 지상국에서는 별도의 원자시계들을 이용하여 시스템 시간을 생성한다. 기존보다 더욱 정확하고 안정된 시스템 시간을 생성하기 위해서는 우수한 성능 (주파수 안정도 및 정확도)을 갖는 세슘원자분수시계나 광시계 같은 새로운 원자시계를 개발하는 것이 필요하다. 그와 동시에 성능은 다소 떨어지지만 이미 상용화 되어 있는 여러 대의 원자시계들을 결합하여 앙상블 (ensemble) 시간을 생성함으로써 어떤 비상 상황에서도 시간을 안정적으로 공급할 수 있도록 시스템의 신뢰도를 높이는 것 또한 중요하다.

앙상블 시간눈금을 가상의 시계 (virtual clock) 또는 문서상의 시계 (paper clock)라고도 부른다. 실물이 있는 시계는 아니지만 컴퓨터 프로그램 또는 데이터 파일(문서)에서 다른 시계들의 기준으로 사용되기 때문이다. 하나의 예로서, 세계협정시인 Universal Coordinated Time (UTC)는 대표적인 가상의 시계인데, 세계 80여개 연구기관에서 생성관리하고 있는 해당 지역의 표준 시간인 UTC(k)와의 시간 차이, 즉 UTC-UTC(k)를 구하는 데 사용된다. 최근에는 차세대 GNSS에 사용할 목적으로 이 가상의 시계를 하드웨어로써 실제로 구현하는 연구가 진행되고 있다 (Godel & Furthner 2017, Godel et al. 2019, Trainotti et al. 2019).

본 논문에서는 앙상블 시간눈금 생성에 널리 사용되고 있는 두 가지 알고리즘, 즉 AT1 알고리즘과 칼만 필터 알고리즘을 비교한 연구 결과를 나타내었다. AT1 알고리즘은 1980년대에 미국 NIST에서 시간눈금 생성 전용 알고리즘으로 개발한 것이다. 현재 국제도량형국 (BIPM)에서 UTC를 생성하는 ALGOS 알고리즘은 AT1의 변형된 형태이다 (Tavella & Thomas 1991, Levine 2012). 이에 비해 칼만 필터는 1960년대 초에 R.E. Kalman이 처음으로 개발한 추정 기법으로 항공, 우주 분야 등에서 널리 사용되고 있는데, 1980년대에 시간눈금 생성에 활용되기 시작했다 (Greenhall 2003). Brown (1991)이 칼만 필터 방법으로 개발한 알고리즘은 현재 Global Positioning System (GPS)에 사용되고 있다.

앙상블 시간눈금을 만들면 앙상블을 구성하는 개별 시계들보다 더 우수한 성능을 얻을 수 있다. 예를 들면, 수소메이저는 단기 안정도는 우수하지만 장기적으로는 표류 (drift) 현상이 나타난다. 이에 비해 세슘원자시계는 단기 안정도는 나쁘지만 장기 안정도는 우수하다. 앙상블 시간눈금은 이런 시계들의 장점을 어느 정도 모두 갖도록 할 수 있다.

앙상블 시간눈금이 이런 특성을 가지려면 개별 시계들의 특성을 잘 반영하도록 적절한 가중치를 부여하는 것이 필요하다. 즉, 안정도가 높은 시계에 더 높은 가중치를 부여하는 것이다. 그런데 안정도는 평균 시간에 따라 달라지므로 시계들의 결합에 특별한 방법이 필요하다. 본 논문에서는 현재 측정 데이터 뿐 아니라 과거 30일간의 데이터를 불러와서 각각 다른 가중치를 부여하는, 이른바 ‘30일 평균 가중치 모델’을 AT1 알고리즘에 사용했다.

한편, Greenhall (2003)은 칼만 필터를 이용한 앙상블 시간눈금 알고리즘을 연구하던 중, 오차 공분산 행렬의 크기가 무한정 증가하는 문제를 해결하는 쉬운 방법을 찾았다. Brown (1991)도 GPS 알고리즘을 개발할 때 이 문제를 이미 알고 있었다. 그래서 그는 공분산의 크기를 주기적으로 줄여주는 방법을 사용했다. 이에 비해 Greenhall은 공분산에서 위상 성분에 해당하는 행렬 값을 0으로 둠으로써 (이것을 x-reduction이라고 부름) 이 문제를 해결했다 (Davis et al. 2005). 이렇게 하더라도 주파수와 주파수 표류에는 아무런 영향을 미치지 않는다는 것을 증명했다. 그는 x-reduction으로 만든 앙상블 시간눈금을 reduced Kalman filter (Kred) 시간눈금이라고 불렀다 (Greenhall 2012). 이 알고리즘에서는 칼만 이득에 내포된 자체 가중치가 자동적으로 각 시계에 부여된다. KRISS에서는 상용의 수소메이저보다 우수한 정확도와 안정도를 갖는 1차 (primary) 주파수 표준기인 세슘원자분수시계 (일명, KRISS-F1, 이하 F1이라 함)를 개발했다. F1를 기준으로 수소메이저의 주파수를 측정하고, 하루 동안의 평균값을 매일 제공한다. 본 논문에서는 두 알고리즘(AT1 및 Kred)에서 생성된 추정 주파수 또는 추정 주파수 표류를 F1으로 보정한 경우와 보정하지 않은 경우를 비교했다. AT1 알고리즘에서는 모든 시계들의 주파수와 주파수 표류를 모두 보정했을 때 가장 좋은 주파수 안정도를 얻었다. 하지만 Kred 알고리즘에서는 기준 시계의 주파수 표류만 보정했을 때 제일 좋은 안정도를 얻었다.

시간눈금 생성에서 사용한 측정 데이터는 온습도가 조절되는 실험실에서 관리하는 실제 시계들의 위상 차이 (시간 차이) 데이터이다. 다시 말하면, 좋은 환경에서 서로 가까이 있는 시계들의 데이터를 이용하기에 측정 잡음을 무시할 수 있는 조건을 갖추었다. 만약 멀리 떨어져 있는 시계들을 원격으로 측정하는 경우에는 잡음을 먼저 필터링하는 과정이 필요하다. 특히 AT1 알고리즘은 측정 잡음을 무시한다는 것을 전제로 구성되어 있다.

알고리즘에 사용된 측정 데이터는 2020년 3월 16일 (MJD 58924)부터 10월 15일 (MJD 59137)까지 약 214일 동안 900초 간격으로 측정한 약 2만개 \(\times\) 3세트의 위상 차이 데이터이다. 이것은 두 대의 수소메이저 (H26, H28)와 두 대의 세슘원자시계 (Cs35, Cs83)에서 H26을 기준으로 한 시간 차이, 즉 H26-H28, H26-Cs35, H26-Cs83를 측정한 것이다.

본 논문의 제 2장에서는 AT1 알고리즘과 Kred 알고리즘을 시계 모델과 연계하여 설명한다. 제 3장에서는 두 알고리즘으로 앙상블 시간눈금 생성 방법을 설명한다. 특히 AT1 알고리즘에서 시계들의 알란분산 그래프로부터 평균 가중치 곡선을 구하고, 가중치를 부여한 방법을 설명한다. 이에 비해 Kred 알고리즘에서는 두 수소메이저의 가중치가 똑 같도록 동일한 잡음 확산 계수를 배정했다. 제 4장에서는 두 앙상블 시간눈금에서 얻은 주파수 안정도 (알란편차 그래프)를 비교, 논의한 후 제 5장에서 결론을 맺는다.

2. TIMESCALE ALGORITHMS

2.1 AT1 Algorithm

어떤 시계가 어떤 시점에서 표시하는 시간은 그 시계가 갖는 특성, 즉 위상 \(x\) (단위: s), 주파수 \(y\) (단위: s/s), 그리고 주파수 표류 \(d\) (단위: s^-1)에 의해 결정된다. 이 세 요소로써 시계의 상태를 나타내는 시계 모델을 3-상태 모델이라고 한다. 그런데 GPS에서는 주파수 표류를 제외한 2-상태 모델을 사용하고 있다. 사용하는 시계의 특성 또는 시계의 모니터링 (또는 조정) 주기 등에 따라 적절한 시계 모델을 사용할 수 있다 (Trainotti et al. 2019).

3-상태 시계 모델에서 시점 \(t\)에서의 시계의 위상 \(x(t)\)는 전 시점 \((t-\tau)\)에서의 \(x, y, d\)를 이용하여 식 (1)과 같이 표현할 수 있다 (Tavella 2008, Levine 2012).

\(x(t)=x(t-\tau)+y(t-\tau)\cdot\tau+0.5\cdot d(t-\tau)\cdot\tau^2  + \xi \)                                                (1)

여기서 \(\xi\)는 랜덤 (random) 잡음을 나타내는데, 평균은 0이고 일정 분산을 가지는 정규분포를 한다고 가정한다. \(\tau\)는 두 시점 사이의 시간 간격이다.

시계의 위상, 주파수, 주파수 표류는 직접 잴 수 있는 양이 아니다. 기준으로 사용한 다른 시계와의 비교에 의해 그 값을 알 수 있다. 그래서 위상 옵셋 (또는 시간 차이), 주파수 옵셋, 주파수 표류 옵셋이라고 부르기도 한다.

AT1 알고리즘에서는 여러 시계들로 만들어진 가상의 앙상블 시계를 기준으로 개별 시계들의 위상 옵셋을 추정하는데, 식 (2)로 표현한다. 여기서 앙상블 시계란 개별 시계들 (앙상블에 포함된 멤버 시계들)의 결합으로 만들어진 시계 또는 시간을 의미한다.

\({\hat x}_{je} (k)=x_{je} (k-1)+y_{je} (k-1)\tau+0.5d_{je} (k-1)\tau^2\)                                             (2)

여기서 \(k\)는 시점\( t_k\)를 나타내며, \(\tau (=t_k-t_{k-1})\)는 측정 시간간격이다. 아래 첨자 \(je\)는 앙상블 시계 \((e)\)에 대한 개별 시계 \((j)\)를 의미하고, 햇 (^) 표시는 추정 값을 나타낸다. 다시 말하면, 위 식은 시점 \(k\)에서 앙상블 시계에 대한 시계 \(j\)의 위상은 이전 시점 \(k-1\)에서의 위상 옵셋, 주파수 옵셋, 주파수 표류 옵셋으로 추정해낸다는 것이다. 이것을 반대로 해석할 수도 있다. 즉 \({\hat x}_{je} (k)\)는 시계 j에 대한 앙상블 시계 \(e\)의 위상 차이를 나타내는 것이다. 이 말은 시계 j를 통해서 앙상블 시간을 추정한다는 뜻이다. 식 (1)이 이론적으로 완벽한 시계를 기준으로 위상 옵셋을 구하는 식이라면 식 (2)는 앙상블 시계를 기준으로 위상 옵셋을 구하는 식이다. 그런데 식 (2)에는 잡음이 반영되어 있지 않다. 측정 잡음이 큰 상황 (예: 원격 측정)에서 측정 데이터를 얻은 경우에는 식 (2)를 적용하기 전에 측정 잡음을 제거하는 과정이 필요하다. 시계 자체가 가지고 있는 잡음은 뒤에서 설명할 시계의 가중치 또는 시 상수 등에 반영하여 처리한다.

시계들 사이의 위상 차이 측정을 위해 앙상블의 멤버 시계들 중에서 하나의 기준 시계를 설정한다. 앙상블에 포함되지 않은 외부 시계를 기준으로 사용할 수 있지만, 시계 수가 충분하지 않은 경우에는 그렇게 하기 어렵다. 앙상블에 포함된 총 시계의 수를 N이라고 하면 측정값은 N-1개가 된다. 식 (2)는 모든 시계에 대해 만들 수 있으므로 식의 수는 N개가 된다. 측정 값이 하나 부족하기 때문에 위 식만으로 앙상블 시간을 바로 알아낼 수 없다.

시점 \(k\)에서 기준 시계 \(r\)과 개별 시계 \(j\) 사이의 위상 차이 측정값을 \(z_{rj} (k)\)로 나타내면, 식 (3)과 같이 쓸 수 있다.

\(z_{rj} (k)=x_r (k)-x_j (k)\)                                                                            (3)

여기서 개별 시계 \(j\)가 기준 시계인 경우에는 \(z_{rr} (k)=0\)이다.

이 알고리즘을 통해 우리가 구하고자 하는 양은 앙상블 시계 \(e\)에 대한 기준 시계 \(r\)의 위상 옵셋이다. 그래서 식 (2)와 식 (3)을 결합하여 식 (4)와 같이 시계 \(j\)에 의한 기준 시계의 위상 옵셋 예측 식을 세운다.

\({\hat X}_{re}^j (k)={\hat x}_{je} (k)+z_{rj} (k)\)                                                                             (4)

모든 개별 시계 \(j\)에 대해 식 (4)를 세울 수 있으므로 식의 수는 N이다.

시계 \(j\)에 부여하는 가중치를 \(w_j\)라고 하면, 식 (4)를 모든 N개 시계에 대해 가중평균하여 식 (5)와 같이 앙상블 시계에 대한 기준 시계의 위상 옵셋 예측 값을 구한다.

\({\hat X}_{re} (k)=\sum\limits_{j=1}^N w_j (k)  {\hat X}_{re}^j (k)\)                                                                           (5)

단, 가중치의 합은 1이다: \(\sum\limits_{j=1}^N w_j (k)=1\).

만약 모든 시계에 가중치가 적절하게 부여되었다고 가정하면 식 (5)의 예측 값은 식 (6)과 같이 실제 값이라고 둘 수 있고, 이 알고리즘의 출력이 된다.

\(x_{re} (k)={\hat X}_{re} (k)\)                                                                                        (6)

여기까지는 NIST에서 사용하는 AT1 알고리즘과 동일하다. 지금부터 본 논문에서 사용하는 방법이 달라지는데, 그 중 가장 큰 차이는 가중치 \(w_j (k)\) 값을 결정하는 방법에 관한 것이다. NIST에서는 다음과 같이 현 시점에서 개별 시계들의 예측 오차 \(\varepsilon_j (k)\)를 구하고, 그 오차의 크기가 과거 일정 기간 동안 오차의 평균 편차에 대해 얼마나 큰지 비교하여 그 상대적인 크기에 따라 가중치를 다르게 부여했다 (Levine 2012).

\(\varepsilon_j (k)={\hat X}_{re}^j (k)-x_{re} (k)\)

이 방법은 매 측정마다 데이터를 검정하고 가중치를 부여한다는 점에서 장점이 있다. 하지만 종류가 다른 원자시계들로 앙상블을 구성한 경우에는 문제가 생긴다. 예를 들면, 수소메이저의 단기 안정도는 상용 세슘원자시계에 비해 수 십 배 더 우수하고, 그래서 예측 오차도 작다. 그 결과, 수소메이저가 대부분의 가중치를 차지하게 된다. 이 때문에 세슘원자시계의 우수한 장기 안정도 특성이 앙상블 시계에 반영되기 어렵다. 그래서 NIST에서는 시계들을 모두 참여시킨다는 목적을 가지고 최대 및 최소 가중치를 정하고, 일정 값 이상 또는 이하가 되지 않도록 가중치를 조정했다. 하지만 본 논문에서는 AT1 알고리즘과 Kred 알고리즘의 성능을 비교하기 위해 가중치에 한계를 두지 않았다.

본 논문에서는 수소메이저와 세슘원자시계의 알란 편차 (ADEV) 그래프를 구하고, 알란 분산 (AVAR)의 평균값에 반비례하는 평균 가중치 그래프를 먼저 구했다. 이 방법은 수소메이저와 세슘원자시계의 서로 다른 두 특성을 효과적으로 결합하는 문제에 봉착했던 미해군관측소 (United States Naval Observatory)가 그들이 개발한 mean timescale 생성에 이용했었다 (Breakiron 1991). 본 논문에서는 수소메이저와 세슘원자시계에 대해 과거 30일 간의 평균 가중치 그래프를 구했고, \(\tau=900\)초 간격의 매 측정마다 과거 30일 데이터 (총 2880개)를 불러와서 과거 시점에 따라 다른 가중치를 부여하여 식 (5)를 계산하고, 식 (6)을 결정했다.

한편, \(z_{rj}\)에 대한 측정이 반복되면 매 사이클에서 식 (2)의 \(x_{je}, y_{je}, d_{je}\)가 업데이트 되어야 한다. 이 중에서 \(d_{je}\)는 \(\tau\)에 대해 일정한 값을 갖는다고 가정한다. 어떤 시계 \(j\)와 앙상블 시계와의 위상 차이 \(x_{je}\)는 식 (6)으로 구한 기준 시계 \(r\)과 앙상블 시계와의 차이 및 측정 값으로부터 식 (7)로 구한다.

\(x_{je} (k)=x_{re} (k)-z_{rj} (k)\)                                                                       (7)

그리고 \(y_{je}\)의 업데이트는 두 단계를 거치는데, 먼저 현 시점에서의 주파수(단위: s/s)는 위상의 차이를 이용해서 식 (8)로 계산한다.

\(f_{je} (k)= {x_{je} (k)-x_{je} (k-1) \over \tau}\)                                                                     (8)

그 후 식 (9)와 같이 exponential filter를 적용하여 과거 주파수는 점차 비율이 작아지도록 한다.

\(y_{je} (k)=y_{je} (k-1)+{\{f_{je} (k)-y_{je} (k-1)\} \over 1+w_y } + d_{je} (k-1) \tau\)                                                 (9)

여기서 \(w_y\)는 무차원 시 상수(time constant)로서 식 (10)으로 결정된다.

\(w_y={T_j \over \tau} \)                                                                                      (10)

여기서 \(T_j\)는 시계 \(j\)의 백색 잡음이 끝나는 평균 시간으로, log-log 알란편차 그래프에서 백색 잡음에 해당하는 \(τ^{-1/2}\)의 기울기가 끝나는 지점의 평균 시간에 해당한다. 일반적으로 \(w_y\) 값은 수소메이저의 경우 수 백, 세슘원자시계의 경우 수 천의 값을 가진다.

AT1 알고리즘은 Fig. 1에 정리된 것처럼 식 (2)부터 식 (10)까지의 과정이 반복되는데, MATLAB을 이용하여 프로그래밍 했다.

Fig. 1. Steps and equations of the recursive AT1 time scale algorithm.

2.2 Kred Algorithm

시간눈금 생성을 위한 칼만 필터 알고리즘은 시간 및 주파수 분야에서 널리 사용되고 있는 알고리즘 중 하나이다 (Greenhall 2003, Galleani & Tavella 2003, 2010). 본 논문에서는 Greenhall의 Kred 알고리즘을 중심으로 기술한다.

Fig. 2. The recursive Kred time scale algorithm.

4대의 원자시계로 구성된 앙상블에 대한 Kred 알고리즘의 단계별 진행은 Fig. 2와 같다. 매 사이클마다 x-reduction을 실시하는데, 식 (11)과 같은 행렬 \(\boldsymbol S\)로 오차 공분산 행렬 \(\Gamma_k\)에 적용한다. 여기서 diag는 블록 대각 행렬 연산자를 의미한다. \(i\)는 시계 번호(1, 2, 3, 4)를 나타내는데, 1번 시계가 기준 시계이다.

\(\boldsymbol S= {\rm diag} (S_1, S_2, S_3, S_4), ~~~~~ S_i=\left [ \begin{matrix}0&0&0\\0&1&0\\0&0&1 \end{matrix}\right]\)                                                       (11)

\(S_i\)의 첫 행 첫 열은 시계 상태 벡터의 위상 항에 해당하는 자리인데, 이것을 0으로 두는 것이다. 오차 공분산 행렬 \(\Gamma_k\)를 \(\boldsymbol {SΓ_k S^T}\)로 대체한다는 것은 모든 시계의 오차 공분산 행렬의 위상 항을 0으로 만든다는 것이다. 이렇게 하더라도 주파수 및 주파수 표류에는 아무런 영향을 미치지 않는다는 것을 Greenhall은 증명했다. 이를 통해 주파수 (위상이 아님)의 분산이 최소화 된다. 그래서 Kred 알고리즘의 출력 중 기준 시계의 추정 주파수 \({\hat y}_1\)가 가장 신뢰할 수 있는 양이다.

Fig. 2에서 첫번째 단계는 시계들 사이의 위상 차이를 측정하는 것이다. 측정값 \(z_{ri}\)는 (기준 시계 \(r\) – 멤버 시계 \(i\))를 의미한다. 그리고 \(k\)번째 사이클에서의 측정값을 \(z_k\)로 나타내었다.

\(k\)번째 사이클에서 앙상블의 상태 벡터는 \(X_k\)로 나타내는데, 이것은 시점 k에서 시계 i의 상태를 나타내는 위상 \((x_i)\) 주파수 \((y_i)\), 주파수 표류 \((d_i)\)로 구성된다. 단, 시계의 상태 벡터는 AT1 알고리즘의 식 (7) 및 식 (9)와 마찬가지로 앙상블을 기준으로 결정된다. 매 사이클에서 구하는 양, 즉 출력은 상태 벡터의 추정 값 \({\hat X}_k\)이다. \({\hat X}_k^-\)는 \({\hat X}_k\)를 구하기 전 예측 값을 나타낸다. 이것은 전 사이클 \((k-1)\)의 상태 벡터로부터 상태 전이 행렬을 적용하여 구한다.

상태 전이 행렬 \(\Phi\)는 개별 시계의 상태 전이 행렬 \(\phi_i\)로 구성되고, \(\phi_i\)는 3-상태 모델일 때 식 (12)와 같다.

\(Φ={\rm diag} (\phi_1, \phi_2, \phi_3, \phi_4), ~~~~~~~ \phi_i=\left [ \begin{matrix} 1&\tau&{\tau^2\over 2} \\ 0&1& \tau \\ 0&0&1 \end{matrix}\right]\)                                                       (12)

오차 공분산 행렬의 예측 값 \(\Gamma_k^-\)도 전 사이클에서 구한 값에 전이 행렬을 적용하여 구한다. 단, 오차 공분산 행렬에는 식 (13-14)와 같이 시스템 잡음의 공분산 행렬 \(\boldsymbol Q\)가 더해진다.

\(\boldsymbol {\Gamma_k^-=\phi\Gamma_{k-1} \phi^T+Q}\)                                                                       (13)

\(\boldsymbol {Q}={\rm diag} (Q_1,Q_2,Q_3,Q_4 ),\)                                                                                                          

\(Q_i=\left [ \begin{matrix} q_1 \tau+{q_2 \tau^3 \over 3} +{q_3 \tau^5 \over 20}&{q_2 \tau^2 \over 2}+{q_3 \tau^4 \over 8}&{q_3 \tau^3 \over 6} \\ {q_2 \tau^2 \over 2}+{q_3 \tau^4 \over 8} & q_2 \tau+{q_3 \tau^3 \over 3} & {q_3 \tau^2 \over 2} \\ {q_3 \tau^3 \over 6} & {q_3 \tau^2 \over 2} & q_3 \tau \end{matrix}\right]\)                                                          (14)

여기서 \(Q_i\)의 구성 요소인 \(q_1,q_2,q_3\)는 각각 White FM, Random Walk FM, Random Run FM 잡음의 확산 계수이다. 시계 종류에 따른 대표적인 값은 (Suess et al. 2010)에 나와 있다.

Fig. 2의 두번째 단계에서 \(\boldsymbol {K_k}\)는 칼만 이득을 나타낸다. \(\boldsymbol {H}\)는 측정 행렬로서 본 논문에서는 4대의 시계 사이의 위상 차이 즉, (기준 시계 – 멤버 시계)를 측정했기 때문에 식 (15) 행렬로 표현된다. 만약 (멤버 시계 – 기준 시계)로 측정한 경우에는 1은 -1로, -1은 1로 바꾸면 된다.

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

칼만 이득 속에 있는 \(\boldsymbol {R}\)은 측정 잡음 공분산 행렬을 나타낸다. Greenhall이 Kred 시간눈금 알고리즘을 처음 만들었을 때는 측정 잡음은 0, 즉 \(\boldsymbol {R}=0\)이라는 가정 하에 만들었다. 그러나 그 후 측정 잡음이 있는 경우에도 Kred 시간눈금 알고리즘이 유효하다는 것을 보였다 (Greenhall 2006). 본 논문에서도 \(\boldsymbol {R}=0\)으로 두고 계산했다.

\(k\)번째 사이클의 최종 상태 벡터 추정 값은 상태 벡터 예측 값 \((\boldsymbol {{\hat X}_k^-})\)에 측정 값과 예측 값의 차이 \((z_k-\boldsymbol {H{\hat X}_k^-})\)에 칼만 이득을 곱한 것을 더해서 구한다: \((\boldsymbol{ {\hat X}_k={\hat X}_k^-+K_k (z_k-H{\hat X}_k^-})\).

이렇게 구한 추정 주파수 및/또는 추정 주파수 표류를 KRISS-F1 주파수 표준기로 보정했다. Kred 알고리즘에서는 기준 시계의 주파수 표류만 보정했을 때의 주파수 안정도가 전체 평균 시간 영역 \((\tau =10^3 \times 10^7 ~{\rm s})\)에서 고르게 우수했다.

3. GENERATION OF TIMESCALES

모의 실험에 사용된 4대의 시계를 이용한 앙상블 시간눈금 생성 시스템의 개략도를 Fig. 3에 나타내었다. 기준 시계 (1번)와 멤버 시계들 (2, 3, 4번) 간의 시간 차이 (위상 차이)를 Dual Mixer Time Difference 방법을 이용하는 Multi-channel Measurement Systems, Microchips (MMS) 장치로써 측정하고, 그 측정값을 데이터 수집 장치 (data acquisition)에 저장한다. 컴퓨터 프로그램에서 이 측정값을 불러와서 알고리즘의 입력 데이터로 사용한다. AT1 또는 Kred 알고리즘으로 계산을 수행하면 출력으로는 앙상블 시계와 개별 시계들 간의 차이, 즉 위상 옵셋, 주파수 옵셋, 주파수 표류의 옵셋이 나온다. 이중에서 위상 옵셋 또는 주파수 옵셋 데이터로 계산한 ADEV 그래프로써 앙상블 시계의 성능 (주파수 안정도)을 평가하고 특성을 분석했다.

Fig. 3. Schematic diagram of the system to generate the ensemble time scale (MMS: Multi-channel Measurement Systems).

Fig. 4. Combined Allan deviation graphs using 3 different data : clock comparison data with respect to the reference clock H26 over the averaging time between 900 s and 86400 s, frequency comparison data with respect the KRISS-F1 between 1 day (86400 s) and 52 day (\(4.5\times10^6\) s), and UTC-clock data from BIPM between 60 day (\(5.2\times10^6\) s) and 120 day (\(1.0\times10^7\) s).

본 논문에 사용된 시계들의 과거 데이터를 이용하여 구한 ADEV 그래프를 Fig. 4에 나타내었다. 이 그래프는 세 가지 다른 종류의 데이터를 이용하여 계산한 것이다. 첫번째는, 수소메이저 H26을 기준 시계로 사용하여 멤버 시계들의 위상 차이 측정 데이터인데, 평균 시간 (\(\tau\)) 900초 \(\sim\) 86400초 (1일) 사이의 그래프에 해당한다. 두번째는, KRISS-F1 주파수 표준기를 기준으로 하루 주기로 측정한 주파수 데이터를 이용하여 계산한 것으로, 평균 시간 86400초 (1일) \(\sim 4.5\times10^6\)초 (52일) 사이의 그래프이다. 마지막 세번째는, BIPM에서 한달 주기로 발표하는 UTC-clock (세계협정시와 시계들의 시간차) 데이터로 계산한 것으로, 평균 시간 \(​ 5.2 \times 10^6 ​\)초 (60일) \(\sim 1.0 \times 10^7\)초 (120일) 사이의 그래프이다. 서로 다른 종류의 데이터를 이용했지만 평균 시간에 대한 ADEV 그래프는 어느 정도 연속성을 보이고 있다. 그림에서 보듯이 단기 안정도는 수소메이저가 세슘원자시계보다 최대 100배 정도 우수하지만 평균 시간 \(\sim 2\times 10^6\)초 이후의 장기 안정도는 세슘원자시계가 더 우수하다.

본 논문에서는 세슘원자시계 Cs35와 Cs83의 그래프를 평균하여 하나의 ADEV 그래프를 얻었다. 그리고 H26에 대한 H28의 ADEV 그래프를 H26과 H28에 공통으로 사용했다.

Fig. 4의 ADEV 값을 제곱하여 얻은 AVAR (알란 분산)에 반비례하는 가중치를 현 시점에서 과거 50일까지 계산한 결과가 Fig. 5이다. 이 곡선은 AVAR 그래프를 fitting하여 구한 3차 다항식으로 계산한 것이다. 각 시점에서 두 종류의 시계에 배당되는 가중치의 합은 1이다. 현 시점에서 수소메이저 (Hm)의 가중치는 최대값인 1에 가깝고, 과거로 갈수록 이것은 점점 줄어들고 있다. 이에 비해 세슘원자시계 (Cs)는 현 시점에서 거의 0이지만 과거로 갈수록 점점 커진다. 과거 30일에 대한 가중치와 측정 데이터를 이용해 계산한 결과와 과거 50일 간의 결과를 비교해 보았는데, 주파수 안정도가 크게 달라지지 않았다. 그래서 가중치 계산에 소요되는 시간 (약 1분)을 고려하여 30일 평균 가중치 모델을 선택했다.

Fig. 5. Weight curves for AT1 algorithm allocated to the hydrogen masers (Hm) and cesium clocks (Cs), which were calculated inversely proportional to the Allan variances (square of ADEV) of Fig. 4.

가중치 외에 AT1 알고리즘에서 결정해야 할 중요한 파라미터는 식 (10)의 시 상수 \(w_y\)이다. 이 값을 결정할 때 과거 실험적으로 구한 ADEV 그래프와 Fig. 4의 ADEV 그래프를 사용했다. \(w_y\) 값을 바꾸어 가며 계산한 각 시계들의 ADEV 그래프 모양과 비교하여 결정했다. 본 논문에서 사용한 시계들의 \(w_y\) 값은 다음과 같다: \(w_{y1}=576,\) \(w_{y2}=384,\) \(w_{y3}=1824,\) \(w_{y4}=5760\).

Kred 알고리즘에서 결정해야 할 중요한 파라미터는 식 (14)에 나와 있는 시스템 잡음의 공분산 행렬의 잡음 확산 계수들 (\(q_1,q_2,q_3\))이다. 그런데 Lee et al. (2020)이 발표한 것처럼 같은 종류의 시계가 3대 이상인 경우에는 3-cornered hat 방법으로 그 시계들의 절대 알란편차 그래프를 구하고 그 모양을 기준으로 fitting하여 확산 계수들을 정할 수 있다. 하지만 현재는 두 대의 수소메이저만이 사용 가능하여 3-cornered hat 방법을 사용할 수 없었다. 그래서 본 논문에서는 두 수소메이저 (H26, H28)의 (\(q_1,q_2,q_3\)) 값을 동일하게 둠으로써 AT1 알고리즘에서 두 수소메이저에 동일한 가중치 그래프를 사용한 경우와 그 결과를 비교해 보기로 했다. 이 값들은 Table 1에 나와 있다. 두 세슘원자시계는 \(q_1\)과 \(q_3\)에서만 미세한 차이가 있다.

Table 1. Dispersion coefficients of the system noise covariance matrix used in the Kred algorithm (AHM: Active Hydrogen Maser, HCC: High-performance Cesium Clock).

Clock (code)	\(q_1\) (unit: s)	\(q_2\) (unit: 1/s)	\(q_3\) (unit: 1/s3)
AHM (H26) AHM (H28) HCC (Cs35) HCC (Cs83)	\(1.0\times10^{-26}\) \(1.0\times10^{-26}\) \(7.0\times10^{-23}\) \(6.0\times10^{-23}\)	\(2.7\times10^{-35}\) \(2.7\times10^{-35}\) \(4.0\times10^{-37}\) \(4.0\times10^{-37}\)	\(4.0\times10^{-51}\) \(4.0\times10^{-51}\) \(3.0\times10^{-53}\) \(4.0\times10^{-53}\)

 

Greenhall은 Kred 시간 눈금을 계산하기에 앞서 Kalman Plus Weight (KPW) 시간 눈금을 먼저 구했다. KPW에서 시계들에 배당되는 가중치는 식 (14)에서 시스템 잡음 공분산 행렬의 1행 1열의 성분, 즉 식 (16)의 \(r_i\) 값에 반비례하여 결정된다 (Greenhall 2003). 여기서 \(i\)는 시계 번호이다.

\(r_i=q_{1i} \tau+{q_{2i} \tau^3 \over 3}+{q_{3i} \tau^5 \over 20}\)                                                                            (16)

이 식으로 구한 KPW의 가중치와 Kred의 칼만 이득에 내포된 가중치는 거의 비슷한 값을 가진다. 그래서 \(12\times3\) 행렬로 나타나는 칼만 이득에 내포된 가중치를 찾는 대신에 위 식으로 가중치를 계산했다.

Table 1의 값과 식 (16)으로 구한 두 수소메이저의 가중치는 각각 0.4996으로 둘의 합은 거의 1이 된다. 그로 인해 세슘원자시계의 가중치는 Cs35의 경우 \(7.1\times 10^{-5}\), Cs83의 경우 \(8.3\times 10^{-5}\)으로 거의 0이다. 가중치가 수소메이저로 편향되는 이 문제는 AT1 알고리즘에서도 예견했던 것이다. 하지만 Kred 알고리즘에서는 칼만 이득의 자체 가중치를 이용했고, Fig. 5의 가중치 그래프를 적용하지 않았다.

KRISS-F1으로 두 알고리즘의 출력 중 주파수 옵셋 \(({\hat y}_i)\) 또는 주파수 표류 옵셋 \(({\hat d}_i)\)을 보정한 방법은 다음과 같다: 하루 동안 (86400초)에 96개의 측정 데이터와 상태 벡터가 생산되는데, 이 중 00:00 UTC (09:00 KST)에 해당하는 주파수 옵셋을 F1의 주파수 값으로 대체한다. 그 외 시간의 95개 주파수 옵셋은 알고리즘의 주파수 출력 값을 사용한다. 그리고 주파수 표류 옵셋의 경우, 00:00 UTC 시점에서 전 날의 주파수와 현재의 주파수의 차이를 86400초로 나눈 값을 사용한다.

F1 데이터로 보정한 결과는 두 알고리즘에서 큰 차이가 났다. AT1에서는 모든 시계의 주파수 및 주파수 표류를 모두 보정해도 주파수 안정도가 크게 좋아지지 않았다. Kred에서는 기준 시계의 주파수 표류만 보정했을 때 장기 안정도에서 크게 개선된 결과를 얻었다.

4. RESULTS AND DISCUSSION

2대의 수소메이저(H26, H28)와 2대의 세슘원자시계 (Cs35, Cs83)로 구성된 앙상블에 Fig. 1의 AT1 알고리즘을 사용하고 Fig. 5의 30일 평균 가중치 모델을 적용하여, 출력으로 앙상블 시계 (TA)에 대한 개별 시계들의 위상 옵셋 \(({\hat x}_{je})\)을 얻었다. 이것으로부터 구한 알란 편차 (ADEV) 그래프는 Fig. 6과 같다. 앙상블 시계에 대한 H26의 ADEV 그래프 (TAH26)와 TAH28의 그래프는 대부분의 평균 시간 영역에서 거의 같은 값을 나타내다가 \(\tau=9.2\times10^5\) s부터 TAH26의 안정도가 좋게 나온다. 이 평균 시간에서부터 TAH26은 세슘원자시계의 우수한 장기 안정도 영향을 받은 것으로 생각된다. \(\tau=2.9\times10^4\) s에서 ADEV=\(7.6\times10^{-16}\) 으로 가장 좋은 안정도를 나타내었다.

Fig. 6. Allan deviation graphs of the 4 clocks, obtained from the phase offsets estimated by the AT1 algorithm.

두 세슘원자시계의 경우, 미세하지만 TACs83이 좀더 좋은 안정도를 나타냈다. 두 시계 모두 동일한 가중치 그래프를 사용했고, 단지 측정 데이터와 시 상수만 다를 뿐이다. 수소메이저의 단기 안정도는 세슘원자시계보다 최대 약 100배 우수하다는 것을 알 수 있다. TAH26의 경우 \(\tau=3.7\times10^6\) s에서의 알란 편차를 제외하면 모든 평균 시간 영역에서 멤버 시계들보다 우수한 안정도를 보이고 있다.

KRISS-F1 주파수 표준기로 모든 시계의 주파수 옵셋 및 주파수 표류 옵셋을 보정했을 때 구한 ADEV 그래프가 Fig. 7에 나와 있다. 두 수소메이저의 장기 안정도가 Fig. 6보다 조금 개선된 것을 알 수 있다. TAH26의 경우, \(\tau=5.8\times10^4\) s에서 ADEV=\(5.4\times10^{-16}\)으로 가장 좋은 안정도를 나타내었다. 특이한 점은 \(\tau=5.8\times10^4 ~ {\rm s} \sim 9.2\times10^5 ~ {\rm s}\) 영역에서는 TAH28이 조금 더 좋은 안정도를 보이고 있다.

Fig. 7. Allan deviation graphs of the 4 clocks, obtained from the phase offsets estimated by the AT1 algorithm. Both the frequency offsets and the frequency drift offsets of the 4 clocks are corrected by KRISS-F1 frequency standard.

AT1 알고리즘은 앙상블 시계를 만들기 위해 개별 시계들을 결합하는 과정이 칼만 필터 방법에 비해 단순하고 직관적이다. 개별 시계들의 성능이 모두 좋아지면 앙상블 시계의 성능이 좋아지기 때문에 F1으로 모든 시계들의 모든 요소를 보정할 때 가장 좋은 안정도를 나타내는 것으로 생각된다.

AT1 알고리즘의 주요 출력은 위상 옵셋이다. 이 알고리즘을 변형시켜 주파수 옵셋을 주요 출력으로 만든 알고리즘이 있는데, AF1 알고리즘이 부른다 (Levine 1996).

Fig. 2의 Kred 알고리즘과 Table 1의 잡음 확산 계수를 적용하여, 그것의 출력으로 얻은 위상 옵셋으로 계산한 ADEV 그래프가 Fig. 8이다. AT1으로 계산한 Fig. 6과 비교해보면 상당히 비슷한 경향을 보인다. 기준 시계 (TAH26)는 \(\tau=5.8\times10^4\) s에서 ADEV=\(8.8\times10^{-16}\)으로 가장 좋은 안정도를 나타내었다. 또한 그 지점부터 TAH28보다 좋은 안정도를 보인다. 두 세슘원자시계도 단기 안정도에서는 Fig. 6과 비슷한데, \(\tau=10^6\) s 이후에는 수소메이저와 비슷한 모양으로 주파수 표류를 보인다는 점이 AT1과 다르다. 이것은 Kred 알고리즘에서는 시계 사이에 상호 작용이 있다는 것을 의미한다. 본 Kred 알고리즘 계산에서 측정 잡음을 0으로 두었기에 측정 값 사이의 상호작용은 무시할 수 있다. 그렇지만 칼만 이득이나 오차 공분산 등의 계산에서 시계들 사이에 서로 영향을 미치는 것으로 생각된다.

Fig. 8. Allan deviation graphs of the 4 clocks, obtained from the phase offsets estimated by the Kred algorithm.

Fig. 9. Allan deviation graphs of the 4 clocks, obtained from the phase offsets estimated by the Kred algorithm. Only the frequency drift offset of the reference clock (H26) is corrected by KRISS-F1 frequency standard.

KRISS-F1 주파수 표준기로 기준 시계 (TAH26)의 주파수 표류 옵셋 \(({\hat d}_1)\)을 보정하여 얻은 ADEV 그래프가 Fig. 9이다. Fig. 8과 달리 TAH26의 장기 주파수 안정도가 대폭 개선된 것을 알 수 있다. \(\tau=5.8\times10^4\) s 이후에서는 대략 \(5\times10^{-16} \sim 9 \times 10^{-16}\)의 범위에서 일정한 주파수 안정도를 나타내고 있다.

그런데 Kred 알고리즘은 위상 옵셋 보다는 주파수 옵셋에 더 최적화된 알고리즘이다 (Davis et al. 2005). Kred에서는 x-reduction을 수행하기 때문에 위상의 절대치보다 위상의 차이 (즉 주파수)가 더 믿을 만한 양이다. 그래서 주파수 옵셋 \(({\hat y}_i)\)과 위상 옵셋 \(({\hat x}_i)\)으로 각각 계산한, 앙상블에 대한 기준 시계 (TAH26)의 ADEV 그래프를 비교한 것이 Fig. 10에 나와 있다. 단, 두 그래프 모두 F1 데이터로 주파수 표류 옵셋을 보정했을 때 구한 것이다. 주파수 옵셋으로 계산한 결과 (yTAH26)는 \(\tau=900\) s에서 ADEV=\(1.5\times10^{-16}\)으로 가장 좋은 안정도를 나타내었다. 이에 비해 위상 옵셋으로 계산한 xTAH26의 경우에는 \(\tau=900\) s에서 ADEV=\(4.2\times10^{-15}\)으로 주파수 옵셋으로 계산한 결과보다 약 28배 나쁘다. 하지만 \(\tau=1.2\times10^5\) s 이후의 장기 안정도는 큰 차이가 없다.

이상에서 살펴본 것처럼 Kred 알고리즘에서는 주파수 옵셋을 이용하는 것이 위상 옵셋을 이용하는 것보다 더 우수한 성능의 앙상블 시계를 만들 수 있다. Fig. 10의 yTAH26 그래프는 앙상블 시계 (TA)와 기준 시계 (H26) 사이의 최고의 주파수 안정도를 나타낸다. 만약 어떤 시계의 주파수를 H26을 매개로 비교 측정한다면 앙상블 시계 (TA)에 대한 그 시계의 주파수 안정도를 계산할 수 있다. 이 작업은 컴퓨터 안에서 이루어진다는 한계를 가지지만 멤버 시계들보다 단기 안정도에서 특별히 우수한 성능을 갖는 앙상블 시계 (가상의 시계)의 유용성을 보여준다.

Fig. 10. Comparison of the Allan deviation graphs obtained from the frequency offset (yTAH26) and the phase offset (xTAH26) estimated by the Kred algorithm. The frequency drift offset of the reference clock (H26) is corrected by the KRISS-F1 frequency standard.

Kred 알고리즘에서 구한 앙상블 시계에 대한 기준 시계의 주파수 옵셋 (\({\hat y}_i\) 또는 yTAH26)은 향후 앙상블 시간을 Micro-Phase Stepper 장치를 이용하여 실제 시계 (주파수 신호)로 구현할 때 활용할 예정이다.

5. CONCLUSIONS

본 논문에서는 4대의 상용 원자시계로 구성된 시계 군에서 개별 시계들보다 우수한 안정도를 갖는 앙상블 시간눈금을 생성하기 위해 세계적으로 널리 사용하고 있는 두 가지 시간눈금 생성 알고리즘, 즉 AT1 알고리즘과 칼만 필터 알고리즘을 비교하는 연구를 수행했다. 온\(\cdot\)습도가 조절되는 실험실에서 관리하고 있는 수소메이저 2대와 세슘원자시계 2대 사이에서 측정한 시간 차이 데이터를 사용하여 두 알고리즘으로 생성한 앙상블 시간눈금의 주파수 안정도를 비교했다. 이 때 시계들에 가중치를 부여하는 방법으로 AT1 알고리즘에서는 30일 평균 가중치 모델을 적용했고, Kred 알고리즘에서는 똑같은 가중치를 갖도록 잡음 분산 계수에 동일한 값을 배정했다. 그리고 KRISS에서 개발한 세슘원자분수시계 (일명, KRISS-F1)를 이용하여 두 알고리즘에서 생성한 주파수 옵셋 또는 주파수 표류 옵셋을 보정한 경우와 하지 않은 경우에 주파수 안정도를 비교했다.

그 결과, AT1 알고리즘에서는 모든 시계들의 주파수 옵셋 및 주파수 표류 옵셋을 보정했을 때 가장 좋은 주파수 안정도가 나왔다. 이에 비해 Kred 알고리즘에서는 기준 시계의 주파수 표류 옵셋 만을 보정한 경우에 가장 좋은 주파수 안정도가 나왔다. 특히 Kred 알고리즘에서 KRISS-F1에 의한 앙상블 시간눈금의 장기 주파수 안정도가 크게 개선되었다. 또한 주파수 안정도 계산에 Kred의 주파수 옵셋 출력을 사용한 경우에는 위상 옵셋 출력을 사용한 경우보다 \(\tau=900\) s에서의 안정도가 약 28배 개선되었다.

AT1 알고리즘에서 30일 평균 가중치 모델을 이용하여 세슘원자시계의 장기 안정도가 앙상블 시간 눈금에 많이 반영되도록 시도했으나 그 효과는 크지 않았다. 또한 Kred 알고리즘에서도 자동적으로 할당되는 가중치의 대부분은 2대의 수소메이저에 한정되었다. 결국 서로 다른 종류의 원자시계보다 한 종류의 원자시계를 사용하는 것이 우수한 앙상블 시간 눈금을 생성하는 데 유리하고, 세슘원자시계보다 수소메이저를 사용하는 것이 더 유리하다는 결론을 얻었다. 그리고 KRISS-F1과 같은 주파수 표준기를 확보하고 있는 경우에는 AT1 알고리즘보다 Kred 알고리즘을 사용하는 것이 장기 안정도 개선에 유리하고, Kred 알고리즘의 출력 중에서 위상 옵셋보다 주파수 옵셋을 사용하는 것이 단기 안정도 개선에 더 유리하다는 것을 알았다.

본 논문을 통해 얻은 이 결과들은 Korea Positioning System (KPS)의 시스템 시간을 생성하는 데 유용하게 사용될 수 있을 것으로 기대한다.