1. INTRODUCTION

정확하고 안정된 시간을 생성하여 보급하는 것은 정보통신 기술의 발전에서 가장 기본적이고 필수적인 요소 중 하나이다. Global Navigation Satellite System (GNSS)를 이용하여 정확한 위치 추정과 항법이 가능한 것도 정확한 기준시간이 공급되기 때문이다. 국제적으로 널리 사용되고 있는 기준시간으로는 '세계협정시', 즉 Coordinated Universal Time (UTC)가 있다. 이 UTC는 세계 여러나라의 80여개 시간관련 연구기관 (또는 시간표준연구실)에서 관리하고 있는 원자시계들을 GNSS위성이나 통신위성을 이용하여 서로 비교 측정한 데이터를 기반으로 생성된다. 이 업무를 총괄하는 국제기구는 프랑스 파리 근교에 위치한 국제도량 형국(BIPM)이다 (Panfilo & Arias 2019).

그런데 UTC는 지연된 시간눈금이다. 다시 말해, 현재 시간을. 알려주는 것이 아니라 첫번째 시각 비교 측정이 이루어진 후 약 45일 뒤에 각 나라의 기준시계가 UTC와 얼마나 차이가 났는지 알려준다. 각 나라의 시간표준연구실에서 유지하는 기준시계 (또는 마스터 시계)의 시간을 UTC(k)라고 부르는데, 이것이 해당 지역 (또는 국가)에서 구현하는 UTC 시간이며 현재 시간을 알려준다. UTC(k)에서 k는 나라 이름 또는 해당 연구기관을 나타내는 4 디짓 이하의 코드로, 우리나라의 경우에는 한국표준과학연구원(KRISS)을 의미하는 'KRIS'이다. 참고로, 현재 운용 중인 모든 GNSS의 시스템 타임은 각 운영 주체국의 국가 표준기관 또는 별도의 목적으로 설립된 국가 기관에서 관장하고 있다. 예를 들면, GPS의 시스템 타임은 미국 해군관측소인 USNO가, Galileo 시스템 타임은 유럽의 여러 표준기관들이 연계하여 생성하고 있다. 이들은 시스템 타임을 UTC에 최상의 안정도(stability)와 불확도(uncertainty)를 갖고 동기 시키는 것을 목표로 한다 (Hahn & Powers 2005) (Lewandowski & Arias 2011).

UTC(k)를 생성하기 위해서는 마스터 시계로 사용하는 원자시계 (일반적으로 세슘원자시계 또는 수소메이저)와 인공위성 시각비교 시스템 및 실험실 내의 시계들을 비교 측정하는 장치 등을 갖추어야 한다 (Lee 2018). 이런 장치를 이용하여 시계들을 5일 간격으로 비교 측정한다. 구체적으로 말하면, Modified Julian Date (MJD)로 4와 9로 끝나는 날의 0시 UTC (한국표준시로 오전 9시 정각)에 시계들을 비교 측정하고, 그 데이터 (즉 시간차이 데이터)를 한달 동안 모아서 인터넷을 통해 BIPM에 제출한다. BIPM에서는 세계 각국의 시간표준연구실에서 보내온 이 데이터를 ALGOS라는 알고리즘으로 분석한 후 매달 Circular T를 통해 [UTC-UTC(k)] 결과를 발표한다. 그래서 해당 월의 첫번째 데이터로부터 약 45일 후에 UTC(k)가 UTC로부터 얼마나 벗어났는지 알 수 있다. 각국의 시간표준연구실에서는 이 시간차이가 일정한 값(±50 ns) 이내에서 유지되도록 마스터 시계를 Micro-phase stepper 등을 이용하여 조정한다.

그런데 [UTC-UTC(k)]가 큰 경우에는 조정해야 하는 시간차이도 커지므로 이를 바로 보정하면 순간적으로 시간변화가 크게 나타난다. 시계를 한꺼번에 많이 조정하면 시계의 안정도가 나빠진다. 그래서 대부분의 시간표준연구실에서는 마스터 시계가 변해온 과거 이력을 바탕으로 미래 어느 시점, 예를 들면, 45일 후에는 어떻게 변했을 것인지 예측하고, 그에 맞추어 매일 (또는 일정 주기로) 조금씩 시간 또는 주파수를 조정한다. 그리고 BIPM 결과가 발표되면 예측했던 값과 실제 값의 차이를 보정하되, 그래도 그 차이가 큰 경우에는 매일 조금씩 나누어 보정한다. 이런 식으로 하면 시간의 흐름에서 갑작스런 변화를 피할 수 있어서 주파수 안정도가 높아진다. 그러므로 안정적인 시간을 지속적으로 유지하기 위해서는 예측을 잘 하는 것이 중요하다. 이를 위해 AT1 알고리즘과 칼만 필터(Kalman filter) 알고리즘이 주로 사용되고 있다 (Levine 2012, Parisi & Panfilo 2016). 이 알고리즘들은 기본적으로 시계 모델을 바탕으로 한다. 즉, 기준시계에 대한 해당 시계의 동기 오차, 주파수 차이 및 주파수 표류 (drift) 값과 이들 각각에 영향을 미치는 잡음을 바탕으로 해당 시계가 미래 어떤 시점에서 어떤 시간을 나타낼 것인지 예측하는 것이다.

이런 알고리즘으로 정확한 예측을 하기 위해서는 시계 자체의 특성을 잘 파악하는 것이 중요하다. 기준시계에 대한 시계의 동기오차나 주파수 차이, 그리고 주파수 표류 같은 양은 사전에 측정하고 분석하면 알 수 있는 양으로, 이것들을 결정론적(deterministic) 양이라고 한다. 이에 비해 잡음은 통계적 (또는 확률적)으로 취급해야 하는데, 시계의 거동 예측을 어렵게 하는 요소다. 그래서 이 논문에서는 잡음 특성을 알아야 하는 시계 모델 대신에 현상학적인 접근법, 즉 일정한 간격으로 측정한 시계의 위상(phase) 데이터만을 이용하여 시계의 거동을 예측하는 방법을 소개한다.

만약 BIPM이 [UTC-UTC(k)] 값을 좀더 자주 발표한다면 각 나라에서는 UTC(k)를 좀더 자주 조정하여 더 정확하게 (UTC에 더 가깝게) 유지할 수 있을 것이다. 여러나라에서 이런 요구가 있었기에 BIPM에서는 '빠른(rapid) UTC'라는 의미인 UTCr 시간눈금을 2013년 7월부터 운영하기 시작했다 (Petit, et al 2014). UTCr은 매일 측정한 시각비교 데이터를 바탕으로 만들어지고, 매주 [UTCr-UTC(k)] 결과가 발표된다. 구체적으로 말하면, 매주 월요일부터 일요일까지 0시 UTC에 비교측정한 데이터를 측정일(D) 또는 늦어도 D+2일 12시 UTC이전에 BIPM의 FTP에 업로드하면 (BIPM FTP 2020), BIPM에서는 첫번째 데이터 측정 10일 후인 수요일 오후(UTC 기준)에 그 결과를 발표한다 (이에 관해서는 Fig. 3에서 자세히 설명함). 2018년 현재 UTCr 생성에 참여하는 연구실 수는 Fig. 1에서 보는 것처럼 62개, 원자시계의 수는 약 300대이다. UTCr은 UTC를 좀더 빠르고 정확하게 구현하기 위해 만들어진 것으로, [UTC-UTCr]의 차이를 줄이는 것을 목표로 한다. 현재 그 시간 변화폭은 Root Mean Square로 약 ±2 ns이다(Panfilo & Arias 2019).

Fig. 1. Time scales generated by BIPM: EAL and TAI are French whose meanings are 'free atomic scale' and 'international atomic time', respectively. TT; Terrestrial Time. PFS/SFS; Primary Frequency Standards/ Secondary Frequency Standards (from Panfilo & Arias 2019).

한편, BIPM에서 사용하는 ALGOS 알고리즘도 업그레이드 되었다. 이 알고리즘은 두 개의 알고리즘으로 구성되는데, 하나는 ‘예측 알고리즘’이고, 다른 하나는 ‘가중치 알고리즘’이다. 예측 알고리즘은 2012년에 업그레이드 되었는데, 원자시계의 주파수 표류를 반영하여 예측할 수 있게 했다 (Panfilo et al. 2012). 그 결과, 세슘원자시계보다 단기 안정도는 훨씬 우수하지만 주파수 표류도 훨씬 큰 수소메이저의 거동을 잘 예측할 수 있게 되었다. 그리고 가중치 알고리즘은 2014년에 업그레이드 되었는데, “예측을 잘 할 수 있는 시계가 좋은 시계다”라는 원칙 하에 수소메이저에 더 큰 가중치를 부여하게 되었다 (Panfilo et al. 2014). 그 결과, BIPM 연차 보고서에 의하면 UTC 생성에서 세슘원자시계는 한 대도 최대 가중치를 받지 못했고, 수소메이저가 최대 가중치를 부여받았다. 이에 따라 UTC(k)를 생성하는 세계 여러 시간표준연구실에서도 마스터 시계로 세슘원자시계 대신에 수소메이저를 이용하는 비율이 늘고 있다. 2018년 현재, UTC 생성에 기여하는 수소메이저의 비율은 약 40%인데, 앞으로 더 늘어날 것으로 예상된다 (BIPM 2018).

본 논문에서는 수소메이저의 거동을 예측하기 위해 일정한 시간간격으로 측정한 데이터, 이른바 시계열(time series) 데이터 분석 방법인 Auto-Regressive Integrated Moving Average (ARIMA) 모델을 이용했다. 본 연구에서는 KRISS에서 운용 중인 2대의 수소메이저가 10일 후에 어떤 값을 나타낼 것인지 예측하고, BIPM에서 매 수요일 오후에 [UTCr-UTC(k)] 결과를 발표하면 예측값과 실제값의 차이를 확인함으로써 예측의 정확도를 확인했다. 이렇게 예측과 확인을 약 100일 동안 반복적으로 실시하여 이 모델의 장점과 단점을 파악했다. 그리고 2대의 수소메이저의 출력 주파수가 서로 강한 상관관계(correlation)를 가진다는 것을 인지한 후 이 상관성을 반영할 수 있는 Vector AutoRegressive Moving Average (VARMA) 모델을 적용하여 두 수소 메이저의 거동을 좀더 정확하게 예측했다. 두 모델을 이용한 분석과 예측을 위해 패키지 프로그래밍 언어인 R을 사용했다. R은 통계 분석에 강점이 있는 여러 패키지들을 무료로 제공하고 있다(The R Project 2020).

UTC 또는 UTCr을 기준으로 UTC(k)를 동기시키는 연구는 각 나라의 시간표준연구실의 주된 임무이다. BIPM 홈페이지(BIPM FTP 2020)에는 UTC 및 UTCr의 생성에 기여하는 세계 각국의 시간표준연구실이 소개되어 있다. 동기시키는 방법에 대해서는 참고문헌 (Whibberley, Davis & Shermar 2011)에 자세히 나와 있다. 하지만 ARIMA 모델을 이용하면 시각동기가 가능하다는 논문은 있으나 (Levine 2012), 저자가 아는 범위에서, 구체적인 방법에 대해서 발표된 것을 아직 찾지 못했다.

본 논문의 제 2장에서는 시계열 자료의 분석 및 예측을 위해 사용한 통계적 모델들인 ARIMA 모델 및 VARMA 모델에 대해 설명한다. 제 3장에서는 시계들의 시간차이 측정 및 데이터 수집 시스템에 대해 설명한다. 제 4장에서는 모델에 따른 예측 결과를 수소메이저와 UTCr의 실제 시간차이 데이터를 이용하여 비교 분석한 후, 제 5장에서 결론을 맺는다.

2. STATISTICAL MODELS

2.1 ARIMA Model

Auto-Regressive Integrated Moving Average (ARIMA) 모델은 1970년대에 Box & Jenkins (1976)가 시계열 분석법들을 통합하여 구축한 것으로, AR 모델, MA 모델, 및 Differencing(차분)과 같은 3가지 작은 모델 (또는 기능)이 결합된 것이다.

AR 모델은 '자기 회귀' 모델이라는 의미로, 현재의 값이 과거의 값에 의존하여 결정되는 경우에 적합하다. 만약 현재 값이 p 시차가 나는 과거의 값에서부터 영향을 받는 경우 이것을 AR(p) 모델이라 부르고 식 (1)로 표현한다 (Shumway & Stoffer 2017).

        \(x_{t}=\phi_{1} x_{t-1}+\phi_{2} x_{t-2}+\cdots+\phi_{p} x_{t-p}+w_{t}\)       (1)

여기서 \(\phi_1, \phi_2, \cdots, \phi_p\) (단, \(\phi_p≠0\))는 AR 모델의 계수이다. w_t는 백색잡음으로 평균은 0이고 분산은 일정하며 정규분포를 갖는다. 그리고 x_t의 평균은 0이라고 가정한다.

MA 모델은 '이동 평균' 모델이라는 의미로, 현재의 값이 과거의 잡음에 의존하여 결정되는 경우에 적합하다. 만약 현재 값이 q 시차가 나는 과거의 잡음에서부터 영향을 받는 경우 이것을 MA(q) 모델이라 부르고 식 (2)로 표현한다.

\(x_{t}=w_{t}+\theta_{1} w_{t-1}+\theta_{2} w_{t-2}+\cdots+\theta_{q} w_{t-q}\)       (2)

여기서 w_t 백색잡음이고, Ø₁, Ø₂, ··· , Ø_q(단, Ø_q≠0)는 MA 모델의 계수이다.

일반적으로 정상성을 갖는 시계열은 AR과 MA 모델을 합친 모델로 설명하는데, 이것을 ARMA(p,q) 모델이라고 부르고, 식 (1)과 (2)를 합한 식 (3)으로 표현한다.

        \(x_{t}=\phi_{1} x_{t-1}+\cdots+\phi_{p} x_{t-p}+w_{t}+\theta_{1} w_{t-1}+\cdots+\theta_{q} w_{t-q}\)       (3)

한편, ARIMA 모델은 비정상(nonstationary) 시계열에도 적용할 수 있다. 이 모델은 차분 기능이 추가된 것인데, 이것이 비정상 시계열을 정상 시계열로 변홖시키는 역할을 한다. 다시 말해, 어떤 시계열을 d차 차분했을 때 ARMA(p,q)로 표현할 수 있다면, 이 시계열은 ARIMA(p, d, q) 모델에 해당하는데, 식 (4)로 나타낼 수 있다.

        \(\phi(B)(1-B)^{d} x_{t}=\theta(B) w_{t}\)       (4)

이 식은 식 (3)에 d차 차분을 합핚 식이다. 식을 간략히 표현하기 위해서 다음과 같은 연산자들을 새로 정의하여 사용핚 것이다.

후향(backshift) 연산자 B는 현 시계열 자료를 한 시차 과거로 변환시키는 역할을 하는 것으로 Bx_t=x_b-1로 표현되고, 이를 중복 적용하면 B²x_t=x_t-2가 된다. 그리고 차분 연산자 ∇는 ∇x_t=x_t-x_t-1 처럼 이웃한 시계열 사이의 차이, 즉 1차 차분을 구하는 연산자이다. 따라서 후향 연산자와 차분 연산자는 ∇x_t =(1-B)x_t의 관계를 가진다. 이것을 두 번 연달아 적용하면, 즉 2차 차분은 식 (5)와 같이 표현된다 (Shumway & Stoffer 2017).

        \(\nabla^{2} x_{t}=(1-B)^{2} x_{t}=\left(1-2 B+B^{2}\right) x_{t}=x_{t}-2 x_{t-1}+x_{t-2}\)       (5)

식 (4)의 좌변에서 (1-B)^dx_t 항은 현 시계열 자료 x_t를 d차 차분하여 정상 시계열이 되도록 만드는 것을 나타낸다. 여기서 연산자 Ø(B)와 θ(B)는 각각 다음과 같이 정의된다.

\(\begin{array}{l} \phi(B)=1-\phi_{1} B-\phi_{2} B^{2}-\cdots-\phi_{p} B^{p} \\ \theta(B)=1+\theta_{1} B+\theta_{2} B^{2}+\cdots+\theta_{q} B^{q} \end{array}\)

2.2 VARMA Model

Vector ARMA (VARMA) 모델은 확률변수가 여러 개인 다변량(multivariate) 시계열의 분석 및 예측에 사용된다. VARMA는 ARMA와 비슷하게 VAR모델과 VMA 모델을 결합한 것이다. 그래서 ARMA와 마찬가지로 정상 시계열에 대해서만 사용할 수 있다. 본 연구에서 계산을 위해 사용한 R 언어의 MTS 패키지에는 비정상 시계열에 적용할 수 있는 VARIMA 함수가 없기 때문에 이와 동일하 기능을 하도록 별도로 Differencing을 수행하여 정상 시계열로 만들고 그 후에 VARMA 모델을 적용했다. 

VAR 모델은 VARMA(p,q) 모델에서 q=0인 경우에 해당하며, 다변량 시계열 분석에서 가장 많이 사용된다. 여기서는 2변량(bivariate)에 국한하여 VAR 모델을 설명한다.

2변량이고 p=1인 경우, 즉 VAAR(1) 모델은 다음식으로 표현된다 (Tsay 2014).

\(z_{t}=\phi_{0}+\phi_{1} z_{t-1}+w_{t}\)

이 식을 벡터와 행렬로 표현하면 식 (6)과 같다.

       \(\left[\begin{array}{l} Z_{a t} \\ z_{b t} \end{array}\right]=\left[\begin{array}{c} \phi_{a 0} \\ \phi_{b 0} \end{array}\right]+\left[\begin{array}{cc} \phi_{a a} & \phi_{a b} \\ \phi_{b a} & \phi_{b b} \end{array}\right]\left[\begin{array}{c} z_{a, t-1} \\ z_{b, t-1} \end{array}\right]+\left[\begin{array}{c} w_{a t} \\ w_{b t} \end{array}\right]\)       (6)

여기서 Ø₀와 w_t는 각각 두 확률변수(z_av, z_bt)의 상수와 백색 잡음을 나타낸다. zat와 zbt는 본 논문에서 각각 시점 t에서 두 수소메이저(Ha, Hb)와 UTCr과 시간차이, 즉 [UTCr-Ha] 및 [UTCr-Hb]를 의미한다. 식 (6)의 행렬에서 비대각(off-diagonal) 요소인 ϕ_ab와 ϕ_ba가 가장 중요한 항이다. 만약 z_at와 z_bt 사이에 상관관계가 없다면 이 두 항은 0이 되고, 위 VAR 모델은 독립적인 2개의 AR 모델로 표현할 수 있다. 구체적인 예를 들면, 만약 이 두 항이 모두 0이라면 두 수소메이저 사이에 교차상관(cross-correlation)이 없다, 즉 두 수소메이저는 독립적이라는 뜻이다.

그런데 ϕ_ab가 0이 아니라면 수소메이저 Ha의 신호는 수소메이저 Hb의 과거 신호 (z_b,t-1)에 의존성을 가진다. 마찬가지로, 0이 아닌 ϕ_ba은 Hb의 신호가 Ha의 과거 신호(z_a,t-1)에 의존하는 것을 나타낸다.

그런데 만약 ϕ_ab=0인데 ϕ_ba≠0이라면, 다시 말해, z_at는 z_b,t-1에 의존하지 않지만 zbt는 z_a,t-1에 의존하는 경우, 영향이 한쪽 방향으로만 미치게 된다. 수소메이저를 예로 들면, Ha는 Hb의 신호에 영향을 미치지만 Hb는 Ha에 영향을 미치지 않는다는 뜻이다.

VMA 모델은 VARMA(p,q) 모델에서 p=0인 경우에 해당하는 모델이다. 2변량이고 q=1인 경우, 즉 VMA(1) 모델은 다음 식으로 표현된다.

\(z_t = \mu + w_t - \theta_1w_{t-1}\)

이 식을 벡터와 행렬로 표현하면 식 (7)과 같다.

       \(\left[\begin{array}{l} z_{a t} \\ z_{b t} \end{array}\right]=\left[\begin{array}{l} \mu_{a} \\ \mu_{b} \end{array}\right]+\left[\begin{array}{l} w_{a t} \\ w_{b t} \end{array}\right]-\left[\begin{array}{ll} \theta_{a a} & \theta_{a b} \\ \theta_{b a} & \theta_{b b} \end{array}\right]\left[\begin{array}{l} W_{a, t-1} \\ w_{b, t-1} \end{array}\right]\)       (7)

여기서 μ는 상수항으로 z_t의 평균을 나타내고, w_t는 백색잡음을 나타낸다. 이 모델에서도 식 (7)의 비대각 요소 θ_{ab​​​​​​}와 θ_ba​​​​​​가 식 (6)에서와 마찬가지로 각각 상대방 확률변수와의 상호작용을 나타낸다.

3. DATA ACQUISITION AND PROCESSING

본 논문에서는 자유동작(free running)하는 두 대의 수소메이저와 UCTr과의 시간차이가 주요 데이터이다. 이 수소메이저들은 BIPM에 등록되어 있는데 BIPM에서 부여한 코드는 각각 405626과 405628이다. 여기서 40은 수소메이저를 의미하고, 나머지 숫자는 등록번호이다. 이 논문에서는 이 등록번호의 맨 뒤 숫자로 구분하여 각각 H6과 H8로 부른다. H6은 2007년 1월에 러시아에서 도입한 것으로 모델명은 Kvarz의 CH1-75A이다. 이것은 2008년 7월부터 KRISS의 마스터 시계로 동작해 왔다. 한편, H8은 2016년 4월에 스위스에서 도입한 것으로, 모델명은 T4Science의 iMaser 3000이다. H6에서 이상이 발생하면서 2019년 8월에 H8로 마스터 시계를 바꾸었다. 본 논문에 사용한 데이터는 H6이 마스터 시계로 동작하던 기간동안 측정된 것이다.

Fig. 2는 UTC(KRIS)를 생성하는 시스템의 개략도이다. 이 그림은 UTC(KRIS)를 위성을 이용하여 외부 기준시와 비교하는 부분과, KRISS 시계 사이의 시간차이를 측정한 후 데이터를 저장하고 분석하는 부분으로 나뉜다. 위성을 이용한 시각비교는 GPS와 같은 항법위성을 이용하거나 통신위성(정지위성)을 이용하여 비교한다. 항법위성은 시스템 타임에 동기된 시각을 전송하는데, 이 시각정보를 이용하여 지상의 기준시계와 비교하고 매 16분마다 데이터를 저장한다(Defraigne & Petit 2015). 통신위성을 이용하는 경우에는 위성에 원자시계가 탑재되어 있지 않으므로 지상의 두 기준시계를 Two-way Satellite Time and Frequency Transfer (TWSTFT) 방법으로 비교한다(ITU 2015). 이렇게 비교한 데이터를 BIPM에 매일 또는 매달 보내면 BIPM은 UTC 생성을 위한 계산에 활용한다.

마스터 시계 H6은 매달 발표되는 [UTC-UTC(KRIS)]의 값에 따라 자체 알고리즘으로 Microphase stepper로써 그 시간차이를 조정(steering)하는데, 그 결과물이 UTC(KRIS)이다. H6 및 H8을 포함한 여러 시계들에서 나온 주파수는 Multi-channel Measurement System (MMS)에 입력되고, 시계들 사이의 시간차이가 출력으로 나와 컴퓨터에 저장된다. MMS에는 UTC(KRIS)의 5 MHz 주파수가 입력되어 각 시계들과의 시간차이를 측정할 때 기준신호로 작용한다.

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

UTCr에 대한 수소메이저 Hn의 시간차이, 즉 [UTCr-Hn]_pred를 예측하는 데 필요한 데이터를 정리하면 다음과 같다.

① UTC(KRIS)-Hn: UTC(KRIS)를 기준으로 Hn (n=6 또는 8)의 시간차이

② UTCr-UTC(KRIS): BIPM에서 매주 발표하는 UTCr과 UTC(KRIS)와의 시간차이

③ [UTCr-Hn]_pred: 본 연구에서 통계 모델로써 예측한 시간차이

①은 매 분마다 측정되어 그 데이터가 컴퓨터에 저장된다. 그중에서 본 논문에서는 매일 UTC 0시에 측정된 값만 사용한다.

②와 ③에 관한 것은 Fig. 3을 참고하여 설명한다. ②는 이번 주(W주) 수요일 밤에, 지난 주, 즉 (W-¹⁾주 월요일부터 일요일까지 7일 간의 데이터가 발표된다. ①과 ②를 더하면, BIPM의 발표가 있는 W주의 수요일 밤에, (W-¹⁾주 월∼일의 7일간에 대한 [UTCr-Hn]_real (실제 측정된 시간차이)을 알 수 있다. 그 때까지의 [UTCr-Hn]_real 데이터를 바탕으로 예측 프로그램을 돌려서 향후 10일간, 즉 (W+¹⁾주의 수요일까지 시간차이를 예측한다. 하지만 예측 당일이 수요일이므로 월, 화, 수 3일에 대한 예측값은 쓸모가 없다. 그러므로 실제로는 W주 목요일부터 (W+¹⁾주의 수요일까지 7일간에 대해 예측한 셈이다. 이렇게 예측한 것이 ③이다. 여기서 예측값 [UTCr-Hn]_pred과 실제값 [UTCr-Hn]_real의 차이를 잔차(residuals)라 부르는데, 이것이 예측 오차(ε)이다. 위의 과정을 매주 반복하면, 수요일 밤마다 지난 주 월∼일에 대한 예측 오차를 확인할 수 있고, 또한 다음 주 목요일까지 예측할 수 있게 된다. 예측 오차, 즉 잔차를 확인함으로써 본 연구에서 사용한 모델들의 성능을 평가할 수 있다.

Fig. 3. Recursive procedure to forecast the (UTCr-clock) data for 7 days from this Thursday to next Wednesday with the [UTCr-UTC(k)] data published by BIPM on this Wednesday night.

만약 Fig. 2에서 나타낸 Microphase stepper를 이용하여 예측값 Δt(≡[UTCr-Hn]_pred)를 0으로 만든다면 그 출력은 UTCr=Hn이 된다. 다시 말하면, Δt 만큼 조정된 Hn이 UTCr과 같다는 뜻이다. 이 출력을 UTC(k)로 사용한다면 UTCr에 근접한 UTC(k)를 얻을 수 있다. 단, ε만큼의 예측 오차가 포함되므로 실제로는 UTC(k) = UTCr + ε이 될 것이다. 결국 UTC(k)를 UTCr에 가깝도록 만들려면 ε를 줄여야 하고, 그러기 위해서는 예측을 잘 해야 한다.

Fig. 4는 MJD 58260(2018.5.22)부터 MJD 58724(2019.8.29)까지 465일 동안의 [UTCr-H6]과 [UTCr-H8]의 변화를 보여준다. 이 기간동안 각각 5201 ns와 4683 ns 만큼 변했다. 이것을 상대주파수 (y=([xi-xi-1]/[ti -ti-1])로 나타내면 각각 1.3×10-13과 1.2×10-13로서, 평균적으로 이런 비율로 H6과 H8의 주파수가 변했다는 의미이다. MJD 58619 (2019.5.16)부터 예측 및 잔차 계산을 시작했다. 다시 말해, 첫번째 10일간의 예측을 위해 그 전 360개의 데이터가 사용되었고, 그 이후부터는 매주 7개씩 추가된 데이터가 예측에 사용되었다. 첫번째 예측을 위해 사용되는 과거 데이터가 충분치 않으면 예측오차가 커지며, 경험적으로 적어도 100개 이상의 데이터가 있어야 데이터 수에 따라 큰 변화가 없는 결과가 나오는 것을 확인했다. 하루에 한개의 데이터를 이용하므로 적어도 100일 이상의 데이터가 축적된 후에 예측을 시작하는 것이 바람직하다.

Fig. 4. Time series data of [UTCr-H6] and [UTCr-H8] for 460 days.

Fig. 4로부터 하루 동안의 상대주파수를 계산하여 그린 것이 Fig. 5이다. 시간차이 데이터에서는 큰 변화가 없어 보이지만 주파수로 그린 것에서는 더 자세한 변화를 볼 수 있다. Y-축의 눈금간격은 1×10-15으로 단위는 무차원(s/s)이다. 그림에서 보는 것처럼 두 그래프는 강한 상관관계를 가진다는 것을 알 수 있다. 이것이 VARMA 모델을 수소메이저의 거동 예측에 사용하게 된 동기이다. 자세한 내용은 제 4장에서 설명한다.

시계열 자료로부터 미래 예측에 적합한 모델을 찾는 방법으로 Box & Jenkins (1976) 방법이 널리 알려져 있다. 이에 대해 간략히 설명하면, 시계열 자료에 대한 자기상관함수(ACF)와 편자기상관함수(PACF)을 구했을 때 시차에 따른 이 함수 값들의 변화 모양으로부터 그에 적합한 모델을 어느 정도 유추할 수 있다. 이 모델의 계수들의 값을 최소제곱법 등으로 추정하고, 그것을 바탕으로 예측값을 구하고 실제값과의 차이(잔차)로부터 모델이 적합한지 검증한다. 예측한 모델이 적합한 경우에는 잔차가 무상관 (uncorrelated) 백색 잡음 형태가 된다. 만약 적합하지 않다면 다시 모델을 바꾸어 앞의 과정을 반복한다. 이러한 반복적인 과정을 효율적으로 처리할 수 있는 프로그램이 요구되는데, 그런 것 중 하나가 R이다. 이 프로그램의 한두 개 함수로써 이 과정들을 한꺼번에 처리할 수 있다. 예를 들면, Fig. 4의 시계열에 적합한 모델을 찾을 때 R에서 'forecast' 패키지를 불러온 후, 함수'auto.arima(x)' (단, x는 시계열 자료)를 사용하면 ARIMA(p,d,q)모델에서 시계열 x에 가장 적합한 p, d, q 의 숫자(정수)와 그 계수들의 값을 찾아준다. 그와 함께 해당 모델의 적합도를 나타내는 잔차의 제곱합(RSS: Residual Sum of Squares)이나 Akaike's Information Criterion (AIC) 값 등도 보여준다.

Fig. 5. Frequency data of [UTCr-H6] and [UTCr-H8], which were calculated from Fig. 4.

본 논문에서 사용한 수소메이저의 경우 p, d, q는 3을 넘는 경우가 거의 없고 대부분 2이내 값을 가졌다. 그 이유로는 수소메이저의 신호가 비교적 단순한 형태로 변한 이유도 있지만 '모수 절약의 원칙'(principle of parsimony)에 따라 가능하면 모수(parameter), 즉 여기서는 p와 q의 값이 작은 것을 선택하도록 함수에 프로그램 되어 있기 때문이다. 이를 인위적으로 해제하는 경우 p와 q 값은 커지지만, 그 결과 (예측 정확도)는 특별한 차이가 없었다. 'Akaike 정보기준'이라는 의미의 AIC는 모델의 적합도를 나타내는 수치인데, 모수가 많을수록 벌점(penalty)을 주어 적 합도가 떨어지도록 정의되어 있다.

적합한 모델을 찾은 후 미래 값을 예측할 때는 forecast 패키지 안에 있는 'forecast' 함수를 사용하되 예측할 데이터 수가 10개라면 'forecast(auto.arima(x), h=10)'을 사용한다.

VARMA 모델의 경우에는 R에서 'MTS' 패키지를 불러오고 그 속에 있는 함수 'VARMA(p,q)'를 사용하면 식 (6)과 (7)에 나타난 것과 같이 상수항, AR 계수 행렬, MA계수 행렬 등을 찾아준다. 단, 이 함수에서는 p와 q 값(정수)을 미리 정해주어야 하는데, 앞에서 설명한 auto.arima(x) 함수의 결과를 이용할 수 있다. 하지만 그 값 주변에서 값을 바꾸어가면서 잔차가 가장 작은 최적 모델(여기서는 계수 값)을 찾을 수 있다. 예측을 위해서는 'VARMApred(h=10)' 함수를 사용한다. VARMA 모델은 정상성을 갖는 데이터만 시계열로 사용할 수 있기 때문에 만약 비정상 시계열이라면 'diff(x)' 함수를 이용하여 차분한 데이터를 사용한다. 1차 차분으로도 정상성을 갖지 못할 경우에는 2차 차분을 실시한다. 차분한 값으로 분석 및 예측을 수행한 후, 원래 데이터 값 (즉, 추세가 포함된 값)으로 예측값을 환원할 때는 'diffinv(x)' 함수를 사용한다. 2차 차분한 경우에는 이 함수를 두 번 적용해야 한다.

Fig. 6. Second difference of [UTCr-H8] of Fig.4 and its auto-correlation function (ACF) and the partial ACF.

Fig. 7. Residuals (=predicted - measured) of [UTCr-H6] and [UTCr-H8]. Other two graphs are shown for comparison.

4. RESULTS AND DISCUSSION

시계열 자료의 정상성(stationarity)은 ARMA 모델이나 VARMA 모델의 전제 조건이다. Fig. 4의 시계열 자료를 1차 차분하면 Fig. 5와 비슷한 모양이 나온다. 차분은 이웃한 데이터 사이의 차이값을 나타낸 것이고, Fig. 5는 그것을 다시 하루(86400초)로 나눈 것이므로 Y-축의 눈금 값과 단위만 달라질 뿐 변하는 모양은 동일하다. 이 그림에서 보는 것처럼 1차 차분한 결과도 기울기를 가지고 있으므로 그 평균과 분산이 시간에 따라 일정하지 않다. 다시 말해, 정상성을 가지지 못하고 있다. 그래서 Fig. 4를 2차 차분했는데 그것을 R의 'tsdisplay(x)' 함수로써 얻은 결과가 Fig. 6이다. 상단 그림에서 기울기가 사라지고 변화폭(분산)도 거의 일정한 것을 볼 수 있다. ACF와 PACF는 시차(lag)에 따라 상관관계가 어떻게 변하는지를 보여준다. 완벽한 백색잡음에 대한 ACF와 PACF는 시차 0에서 1이 나오고, 0이 아닌 시차에서는 모두 0이 나온다. 하지만 실제 상황에서 완벽한 백색잡음은 존재하지 않기 때문에 일정한 범주(예: 그림의 점선 또는 5%) 내에 들면 정상성을 가진 것으로 본다. ACF와 PACF의 모양으로부터 이 시계열에 적합한 모델을 추정할 수 있다. 이 시계열 자료에 대해 auto.arima(x) 함수를 적용한 결과 최적 모델은 ARIMA(p=0,d=2,q=2)가 나왔다. 다시 말하면, Fig. 4를 2차 차분한 시계열 자료는 MA(2) 모델에 가장 적합하다는 것이다.

105일 동안 15 차례에 걸쳐 반복적으로 두 수소메이저 (H6, H8)와 UTCr의 시간차이를 예측하고, 실제값과의 차이, 즉 잔차를 그린 것이 Fig. 7의 가운데 부분에 있는 두 그래프이다. 이 기간 동안 최대 잔차폭 (최고값-최저값)는 각각 6.1 ns와 6.2 ns가 나왔다. 이 두 그래프는 제 3장에서 설명한 것처럼 ARIMA 모델로써 예측한 [UTCr-Hn]≡Δt 값을 Fig. 2의 Microphase stepper로써 상쇄시키는 경우 남는 예측오차 ε의 크기를 나타낸다. 단, 이 때 Microphase stepper의 기능 때문에 시간차이(즉 위상)를 바로 보상하는 것이 아니라 5 MHz 주파수를 변화시켜 보상한다. 이 결과와 비교하기 위해 BIPM에서 발표한 UTC-UTC(KRIS)와 UTCrUTC(KRIS)를 같이 그렸다. 이 기간동안 최대 변화폭은 약 37 ns로서, 예측 오차의 변화폭에 비해 최대 6배 크다.

Fig. 4의 두 그래프, 즉 [UTCr-H6]≡T6과 [UTCr-H8]≡T8 사이의 상관관계를 알아보기 위해 T6과 T8의 1차 차분 및 2차 차분에 대한 교차 상관(cross-correlation)을 구해 보았다. 식 (5)에서 사용한 차분 연산자 ∇를 사용하면 각각 다음 식으로 표현된다. 단, 두 식에서 n은 수소메이저를 나타내는 것으로 n=6 또는 8이다.

\(\nabla T n=[U T C r-H n]_{t}-[U T C r-H n]_{t-1}\)

\(\nabla^{2} T n=(\nabla T n)_{t}-(\nabla T n)_{t-1}\)

∇T6과 ∇T8 사이의 교차 상관을 R의 함수 'ccf(x)'를 이용하여 구한 결과가 Fig. 8a이다. 1차 차분을 했지만 여전히 강한 상관관계를 가진다는 것을 알 수 있다. 2차 차분한 ∇² T6과 ∇² T8 사이의 교차 상관을 구한 결과가 Fig. 8b이다. 아래쪽 네모 박스로 표시한 시차 +1과 -1을 제외하면 다른 시차에서는 상관관계가 대폭 줄어든 것을 알 수 있다. 현 시점과 가장 가까운 ±1 시차에서는 여전히 상관관계가 남아 있다는 것을 뜻한다.

Fig. 8. Cross-correlation of [UTCr-H6] and [UTCr-H8]: (a) their first differences and (b) their second differences.

2차 차분한 시계열 자료를 이용하여 VARMA 모델로써 예측하였고, 원자료에 대한 ARIMA 모델의 결과와 비교했다. MJD 58672부터 58713까지 42일(=6주)에 대해 예측과 잔차 확인을 반복한 결과, [UTCr-H6]의 잔차는 Fig. 9와 같다. ARIMA(p=0,d=2,q=2)와 diff=2 & VARMA(p=0,q=2)의 최대 잔차를 비교해보니 양쪽 모두 2.9 ns가 나왔다.

Fig. 9. Comparison of results by ARIMA and VARMA models for [UTCr-H6].

마찬가지 방법으로, [UTCr-H8]의 잔차를 구한 결과는 Fig. 10과 같다. 최대 잔차는 ARIMA의 경우는 4.5 ns, VARMA의 경우는 2.0 ns로 VARMA 모델이 더 좋은 예측 결과를 보였다.

Fig. 10. Comparison of results by ARIMA and VARMA models for [UTCr-H8].

H6의 경우에는 VARMA나 ARIMA가 같은 결과를 보이지만 H8의 경우에는 VARMA가 ARIMA보다 더 좋은 결과를 보였다. 이것은 식 (6)과 (7)에서 설명한 상관관계에 따른 상호작용에서 H8은 H6의 좋은 영향을 많이 받았지만 H6은 H8로부터 별 영향을 받지 않았다는 것을 알 수 있다. 이에 대해서는 뒤에서 좀더 자세히 살펴본다.

2차 차분한 시계열에 대해 VARMA(p=0, q=2)가 최상의 결과(제일 작은 잔차)를 생산한 것인지 여부를 확인하기 위해 p와 q값을 바꾸면서 잔차의 변화폭을 비교해봤다. VARMA(p=1, q=1)에 의한 [UTCr-H6]의 최대 잔차폭은 2.8 ns였고, [UTCr-H8]은 2.0 ns로서, 앞의 2.9 및 2.0 ns와 거의 같은 결과를 보였다. VARMA(p=1, q=2)에서는 각각 2.9 ns와 2.4 ns가 나왔다. 결론적으로 auto.arima 함수가 제시하는 p=0, q=2를 사용해도 좋다고 볼 수 있다.

본 논문에서처럼 두 수소메이저가 동일한 잡음 형태 (즉 p, q 값)를 가지는 경우에는 VARMA 모델을 사용하는 데 별 어려움이 없다. 하지만 두 수소메이저가 전혀 다른 값, 다시 말하면 서로 다른 잡음 형태를 가지는 경우에는 최적의 p, q 값을 정하는 데 어려움이 있다. 그리고 이런 경우에는 둘 사이에 상관성도 높지 않을 것으로 예상되고, VARMA를 이용한 장점이 나타나지 않을 수도 있다. 본 논문의 두 수소메이저는 뒤에서 설명한 것처럼 강한 상관관계를 갖고 있는데, 그로 인해 동일한 p, q 값을 가지는 것으로 여겨진다.

다변량 변수의 해석에서 충격반응함수(Impulse Response Function)가 있다. 이는 서로 상관관계를 가지는 변수들 사이에서 어느 쪽이 더 큰 영향을 미치는지 알아볼 때 쓰는 함수이다. 예를 들면, 동일한 시기에 측정한 여러 나라의 경제 지표(예: GDP)가 있을 때 어느 나라가 어느 나라로부터 더 큰 영향을 받는 지 등을 분석하는 데도 사용된다(Tsay 2014).

본 논문에서는 R의 'VARMAirf' 함수를 이용하여 H6과 H8 사이에서 충격반응을 살펴보았다. 충격에 대한 반응이란, 예를 들면, VARMA 모델, 즉 식 (6)과 (7)이 합쳐진 식에서 zat (또는 zbt)를 t=0 시점에서 1만큼 값을 증가시켰을 때 zt+j(단, j>0)에서 어떤 변화가 나타나는 지 본다는 뜻이다. 그 결과 Fig. 11의 그래프들을 얻었다. 그림의 X-축은 미래(t>0)시차를 나타내고 Y-축은 IRF 함수 값을 나타낸다. Y 값이 높다는 것은 충격의 반응이 크다는 것을 뜻한다. 상단 오른쪽 그림은 H8이 H6에 미치는 영향을 나타내는데 전반적으로 값이 작고 시차가 커질수록 (미래로 갈수록) 줄어들고 있다. 이에 비해 하단 왼쪽은 H6이 H8이 미치는 충격을 나타내는 것으로 시차에 따라 점점 증가하는 것을 알 수 있다. 반면에 하단 오른쪽은 H8이 자기 자신에게 미치는 영향을 나타내는데 시차가 커질수록 급격히 줄어든다. 이 결과는 앞서 VARMA와 ARIMA의 결과로부터 추정한 것과 동일하다. 다시 말해, H6이 H8에 강한 영향을 미치면서 H8의 예측값이 더 좋아진 것이다. 이것은 곧 H6이 더 안정적으로 동작하고 있다는 것을 의미하는데, Fig. 4와 Fig. 5의 그래프에서도 이런 결과를 어느 정도 추정할 수 있다. 하지만 왜 H8이 H6에 나쁜 영향은 미치지 않고, H6이 H8에 좋은 영향만 미치는 지에 대해서는 정확히 알 지 못한다. 단지 다음과 같은 점이 원인일 수도 있을 것이다: H6의 동작상태를 나타내는 지표 값이 한동안 정상범위를 벗어나서 MJD 58724(2019.8.29)에는 수리를 위해 마스터 시계를 교체하게 되었다. 이로 인해 마스터 시계로 동작하는 기간 동안에 잡음의 세기가 H8보다 더 강했기 때문일 수 있다.

Fig. 11. Results of the impulse response function (IRF) between H6 and H8. The arrow indicates the direction of influence

두 수소메이저가 상관성을 갖게 된 원인에 대해서는 측정시스템, 구체적으로 Fig. 2의 MMS가 원인일 수 있다는 논문이 있다(Levine 2012). 시간차이를 정밀하게 측정하는 방법으로 'dualmixer time difference' 방법이 있다. 이것은 측정 장치 속에 별도의 발진자를 사용하여 측정하고자 하는 두 주파수원 (여기서는 H6와 H8)과 로컬 발진기 (local oscillator) 사이의 차이 주파수, 즉 더 낮은 주파수 2개를 만듦으로써 두 주파수 차이, 곧 시간차이를 더 정밀하게 알아내는 방법이다. 이 때 H6과 H8의 주파수가 로컬 발진기와 동시에 비교될 때 두 신호가 상관관계를 갖게 된다는 것이다. MMS도 이런 방식으로 동작하기 때문에 두 수소메이저가 상관관계를 가지게 된 것으로 추정된다.

5. CONCLUSIONS

본 논문에서는 [UTC(KRIS) + 9 시간]으로 정의된 한국표준시를 생성하는데 사용되고 있는 두 수소메이저(H6과 H8)와 UTCr과의 시간차이, 즉 [UTCr-H6]과 [UTCr-H8]의 미래 값을 예측함으로써 더 정확한 UTC(KRIS)를 만드는 방법에 대해 논하였다. 예측을 위해 통계 모델인 ARIMA 와 VARMA 모델을 적용했으며, 패키지 프로그래밍 언어인 R을 사용하여 계산했다. ARIMA 모델을 적용했을 때 예측 정확도는 잔차의 최대폭으로 나타낼 때 약 6.2 ns로서, UTC-UTC(KRIS)의 최대 변화폭에 비해 약 6분의 1 수준이었다. 다른 말로 하면, 6배 더 정확한 UTC(KRIS)를 생성할 수 있는 가능성이 있다.

H6과 H8 사이에 강한 상관관계를 고려하여 VARMA 모델을 사용하여 예측한 결과, [UTCr-H8]의 예측값은 잔차의 최대폭이 ARIMA 모델 (4.5 ns)에 비해 절반 이하(2.0 ns)로 줄어들었다. 하지만 [UTCr-H6]은 ARIMA와 비슷한 결과(2.9 ns)가 나왔다. 충격반응함수(IRF)를 적용해본 결과, H8은 H6으로부터 강한 영향을 받지만 H6은 H8의 영향을 거의 받지 않는다는 것을 확인할 수 있었다. 두 시계가 상관관계가 있고 또 잡음 형태가 비슷한 경우에 VARMA 모델로써 더 정확한 예측이 가능하다는 것을 알 수 있었다.

결론적으로, ARIMA 모델을 적용하여 UTC(KRIS)을 생성하게 되면 정확도를 현재보다 6배 더 높일 수 있고, VARMA 모델을 적용하면 그 보다 더 개선할 수 있음을 알 수 있었다. 그리고 본 논문에서 구해진 결과는 향후 구축 예정인 한국형 위성항법시스템(KPS: Korean Positioning System)과 같은 고정밀의 시스템 타임이 요구되는 항법시스템에 활용될 수 있을 것으로 기대된다.

AUTHOR CONTRIBUTIONS

Methodology and software, H.S. Lee; data acquisition, T.Y. Kwon; measurement system, Y.K. Lee, S.-h. Yang, D.-H. Yu.

CONFLICTS OF INTEREST

The authors declare no conflict of interest.