Carrier Phase Based Cycle Slip Detection and Identification Algorithm for the Integrity Monitoring of Reference Stations

Abstract

In order to ensure the high-integrity of reference stations of satellite navigation system, cycle slip should be precisely monitored and compensated. In this paper, we proposed a cycle slip algorithm for the integrity monitoring of the reference stations. Unlike the legacy method using the Melbourne-Wübbena (MW) combination and ionosphere combination, the proposed algorithm is based on ionosphere combination only, which uses high precision carrier phase observations without pseudorange observations. Two independent and complementary ionosphere combinations, Ionospheric Negative (IN) and Ionospheric Positive (IP), were adopted to avoid insensitive cycle slip pairs. In addition, a second-order time difference was applied to the IN and IP combinations to minimize the influence of ionospheric and tropospheric delay even under severe atmosphere conditions. Then, the cycle slip was detected by the thresholds determined based on error propagation rules, and the cycle slip was identified through weighted least square method. The performance of the proposed cycle slip algorithm was validated with the 1 Hz dual-frequency carrier phase data collected under the difference levels of ionospheric activities. For this experiment, 15 insensitive cycle slip pairs were intentionally inserted into the raw carrier phase observations, which is difficult to be detected with the traditional cycle slip approach. The results indicate that the proposed approach can successfully detect and compensate all of the inserted cycle slip pairs regardless of ionospheric activity. As a consequence, the proposed cycle slip algorithm is confirmed to be suitable for the reference station where real time high-integrity monitoring is crucial.

1. INTRODUCTION

4차 산업혁명이 도래함에 따라 사용자에게 고정밀, 고신뢰성 위치(Positioning) 및 항법(Navigation), 시각(Timing) 정보를 제공해주는 위성항법 시스템은 미래 국가 성장 동력의 핵심으로 부각되고 있다. 이에 따라 우리나라에서도 한반도 주변 사용자를 위한 지역 위성항법 시스템인 한국형 위성항법 시스템(Korean Positioning System, KPS)을 개발 중이다. 위성항법 시스템의 중요 핵심 알고리즘인 위성 시계, 궤도 오차 추정, 그리고 전리층 지연 오차 추정을 위해서는 기준국 데이터의 무결성이 보장이 중요하다. 무결성이 보장되지 않는 데이터를 기반으로 생성된 궤도력과 전리층 모델의 사용은 위성의 Signal In Space - User Range Error에 영향 주어 결과적으로 위성항법 시스템의 성능을 전반적으로 저하 시키게 된다. 따라서 무결성이 보장된 기준국 데이터만을 사용하도록 하는 무결성 감시 과정이 필수적이다. 위성항법 시스템의 무결성 감시는 크게 이상 검출과 오차 추정의 두가지 기능으로 나눌 수 있으며, 이중 이상 검출은 시스템 내에서 발생할 수 있는 이상을 검출하여 이를 제거하거나 사용하지 않도록 하기 위한 기능으로 Signal Quality Monitoring, Data Quality Monitoring (DQM), Measurement Quality Monitoring 등이 해당된다 (Yun 2007). 이러한 이상 검출 알고리즘들은 기본적으로 모니터링 값을 계산하고, 주어진 무결성 요구 조건을 기반으로 계산된 임계값과 비교하는 방법으로 측정치의 이상을 판단하고 제거하는데, 이를 위해서는 반송파 데이터에 포함된 사이클 슬립(Cycle-Slip)의 검출 및 보상이 전처리 과정에서 필수적으로 선행되어야 한다 (Bisnath et al. 2001). 사이클 슬립은 위성과 수신기 사이에 존재하는 사이클의 초기 위상정수를 의미하는 모호정수에 발생하는 순간적인 점프 현상을 의미한다. 모호정수는 위성신호 수신에 단절이 발생하지 않는 한 초기에 설정된 일정한 미지정수로 유지되는데, 신호의 단절이 발생하는 경우, 수신기의 위상정수 계산기가 초기화 되어 정수의 사이클의 수에 의해 누적되는 모호정수에 순간적인 점프 현상이 발생한다. 사이클 슬립의 발생은 위성에서부터 수신기 사이의 거리를 예측 불가능하게 하여 기준국 데이터의 무결성에 심각한 영향을 미치므로 연속적이고 정확한 사이클 슬립의 감시 및 식별이 요구된다.

사이클 슬립 검출 및 보상에 관한 다양한 연구가 과거부터 진행되어오고 있다. 특히, 이중 주파수를 이용한 사이클 슬립 검출 알고리즘이 많이 제안되었다. Goad (1985)는 Geometry-Free (GF) 선형조합을 통해 이중 주파수 측정치 내에 포함된 공통 오차 성분을 제거하고, 전리층 오차의 변화율을 이용하여 사이클 슬립을 탐지하는 Phase Ionospheric Residual (PIR) 알고리즘을 처음 제안하였다. PIR과 유사한 방법으로 Liu (2011)는 이중 주파수 관측값의 GF 선형조합을 이용하여 계산된 Total Electron Content 변화율을 이용하여 사이클 슬립을 검출하였다. 또한, Cai et al. (2013)은 PIR에 포함된 전리층 영향을 최소화하기 위하여 2차 시간 차분을 추가로 적용한 Second-order, Time-difference PIR을 제안하였다. Blewitt (1990)은 GF 선형조합과 더불어 Melbourne-Wübbena (MW) 선형조합 (Melbourne 1985, Wübbena 1985)을 사용하여 사이클 슬립을 검출하는 TurboEdit 알고리즘을 처음 소개하였다. MW 선형조합을 사용하면 간단히 모호정수 항을 분리할 수 있으며, Wide-lane (GPS L1, L2의 경우 86.2 cm)이 생성되어 모호정수 해석에 효과적인 선형조합이다. 또한, MW 선형조합을 사용하면 GF 선형조합에 의해 탐지가 불가능한 일부 사이클 슬립 조합의 탐지가 가능하므로 GF 선형조합과 MW 선형조합을 함께 이용하면 보다 효과적인 사이클 슬립 탐지가 가능하다. 그러나 반송파 관측값 만을 사용하는 GF 선형조합과는 다르게 MW 선형조합은 정확도가 낮은 의사거리 관측값이 추가적으로 사용된다. 따라서 MW의 경우는 작은 크기의 사이클 슬립이 발생하였을 때, 탐지가 쉽지 않다는 한계가 있다. 따라서 Cai et al. (2013)은 MW 선형조합에 포함된 관측 오차의 영향을 완화하기 위하여 전후방 스무딩 필터를 적용하는 Forward and Backward Moving Window Averaging 알고리즘을 제안하였다. 그러나, 해당 알고리즘은 후처리용으로 만 사용 가능하므로 실시간 사이클 슬립 감시가 요구되는 경우에는 적용이 어렵다는 문제가 있다. Kim et al. (2018)은 Real-Time Kinematic (RTK) 기준국 데이터의 무결성 감지를 위하여 전리층 기반 선형조합을 사용하는 새로운 사이클 슬립 검출 및 보상 알고리즘을 제안하였다. 해당 논문에서 저자는 상호 보완적으로 사용할 수 있는 반송파 데이터 기반 IN, IP 선형조합을 적용하여 사이클 슬립 탐지 성공율을 높힐 수 있었다. 또한, 이차 시간 차분을 적용하여 전리층과 대류층 지연에 의한 영향이 최소화된 모니터링 값을 분석하여 오탐지율을 크게 낮추었다.

본 논문에서는 기준국 데이터의 무결성 검사 중 하나인 DQM에서 중요한 기준 값으로 사용되는 반송파 측정치의 무결성 감시를 위한 고성능 사이클 슬립 검출 및 보상 알고리즘 제안한다. 기준국 수신기의 정확한 위치 정보 및 방송 궤도력 등의 정보를 기반으로 반송파 내에 포함된 각종 오차 요소를 제거하는 전처리 과정이 우선 수행된다. 이어서 보정된 반송파 데이터의 IN, IP 선형조합을 기반으로 사이클 슬립 검출 과정을 수행하며, 이때 전리층과 대류층 오차의 영향을 최소화하기 위하여 2차 시간 차분된 값을 모니터링 값으로 사용한다. 또한, 각 주파수 반송파 데이터에 포함된 사이클 슬립의 크기를 결정하기 위해서는 IN, IP 선형조합의 모니터링 값을 이용한 가중 최소 제곱법 (weighted least squares)을 적용한다. 이러한 일련의 과정을 통해 높은 정확도의 반송파 데이터만을 이용한 고성능 사이클 슬립 검출 및 보상이 가능하다. 본 논문에서 제안하는 사이클 슬립 알고리즘의 성능 평가를 위해서 전리층 비활동기 및 활동기 데이터를 이용한 실험이 진행되었다.

2. CARRIER PHASE RESIDUAL

본 논문에서 제안하는 사이클 슬립 알고리즘은 이중 주파수 수신기로부터 수신된 GPS L1, L2 반송파 관측값을 기본으로 사용한다. 위성 \(i\)로부터 수신된 L1, L2 반송파 관측값은 Eq. (1)과 같이 표현된다.

\(Φ_1^i=ρ^i+c(dt-dT^i )-I+T+λ_1 N_1^i+ε_{Φ_1^i}\)                                                                  
\(Φ_2^i=ρ^i+c(dt-dT^i )-γI+T+λ_2 N_2^i+ε_{Φ_2^i}\)                                                           (1)

여기서 \(Φ\)는 수신기에서 관측된 반송파 관측값이며, \(ρ\)는 위성과 수신기 사이의 기하학적 거리이다. \(dt\)와 \(dT\)는 각각 수신기와 위성의 시계오차를 나타내며, \(c\)는 빛의 속도(299792458 m/s\(^2\))이다. 또한, \(I\)는 위성에서 출발한 L1 신호가 전리층을 통과하면서 발생하는 신호의 지연을 나타내는데, L1 (1575.42 MHz), L2 (1227.6 MHz) 주파수 비율의 제곱인 \(γ(f_1^2/f_2^2)\)를 곱하면 L2 신호의 전리층 지연값을 계산할 수 있다. \(T\)는 대류층에서 발생하는 신호의 지연을, \(N\)은 모호정수를, \(λ\)는 각 주파수의 파장을 의미한다.

본 논문은 기준국의 무결성 감시를 위한 사이클 슬립 탐지 알고리즘을 제안한다. 즉, 사용자는 이미 기준국 수신기의 정확한 위치를 알고 있으며, 위성에서 방송하는 궤도력을 이용하면 위성의 위치 및 위성의 시계오차의 계산이 가능하다. 수신기와 위성의 위치를 이용하여 계산된 기하학적 거리(\(ρ\))와 위성의 시계오차(\(dT\))를 반송파 관측값에서 보정하면 Eq. (2)와 같이 표현할 수 있다.

\(\tilde{Φ}_1^i=λ_1 Φ_1^i-ρ^i+cdT^i=cdt-I+T+λ_1 N_1^i+ε_{\hat{Φ}_1^i}\)                                                      
\(\tilde{Φ}_2^i=λ_2 Φ_2^i-ρ^i+cdT^i=cdt-γI+T+λ_2 N_2^i+ε_{\hat{Φ}_2^i}\)                                              (2)

\(Δ\tilde{Φ}_1^i (t)=\tilde{Φ}_1^i (t)-\tilde{Φ}_1^i (t-1)=cΔdt-ΔI+ΔT+λ_1 ΔN_1^i+ε_{Δ\hat{Φ}_1^i (t)}\)                                    
\(Δ\tilde{Φ}_2^i (t)=\tilde{Φ}_2^i (t)-\tilde{Φ}_2^i (t-1)=cΔdt-γΔI+ΔT+λ_2 ΔN_2^i+ε_{Δ\hat{Φ}_2^i (t)}\)                             (3)

Eq. (3)은 보정된 반송파 관측값의 시간 차분 결과를 보여준다. 여기에서 모호정수의 변화량(\(ΔN\))은 사이클 슬립의 크기를 의미하며, \(Δdt\)는 수신기의 시계 drift를 나타낸다. 또한, \(ΔI\)와 \(ΔT\)는 각각 전리층과 대류층 지연의 변화 속도를 나타내는데, 일반적으로 전리층과 대류층 지연은 시간에 따라 매우 천천히 변화한다고 알려져 있다 (Park et al. 2006, Olynik et al. 2002, Olynik 2002). Park et al. (2006)의 International GPS Service (IGS) 대류층과 전리층의 천정방향 지연값을 분석한 결과에 의하면 대류층 지연은 최대 2.2 mm/s의 속도로 변화하며, 전리층 지연의 경우는 최고 2.03 cm/s의 속도로 변화한다. 또한, Olynik et al. (2002)과 Olynik (2002)이 이중 주파수 관측값을 이용하여 계산한 전리층 지연값을 분석한 결과에 따르면 전리층 폭풍이 발생하지 않은 환경에서 전리층 지연의 변화율은 0.6에서 1.2 cm/s로 Park et al. (2006)과 유사한 변화율을 나타내었다. 즉, 전리층 폭풍이 발생하지 않은 일반적인 대기 상태에서 Eq. (3)의 시간 차분 관측값에 포함된 전리층과 대류층 지연의 변화율은 사이클 슬립(GPS L1 파장: 19.0 cm, L2 파장: 24.4 cm)과 비교하여 매우 작은 값이라고 할 수 있다.

\(c\cdot Δ\hat{dt}={1\over m} \sum\limits_{i=1}^m Δ\tilde{Φ}_{IF}^i\)                                                                                             (4)

\(Δ\tilde{Φ}_{IF}^i=a_1 Δ\tilde{Φ}_1^i-a_2 Δ\tilde{Φ}_2^i=CS_{IF}+c\cdot Δdt+ΔT+ε_{Δ\hat{Φ}_{IF}^i}\)                                               (5)

수신기의 시계 drift의 경우는 Eqs. (4-5)와 같이 Ionosphere-Free (IF) 선형조합의 평균으로 계산 가능하다. Eq. (5)의 \(a_1\)과 \(a_2\)에 각각 \(γ\over γ-1\)과 \(1\over γ-1\)을 적용하여 생성되는 IF 선형조합을 이용하면 전리층 항의 완벽한 소거가 가능하며, 대류권 지연의 변화 속도\((ΔT)\)의 경우는 작은 값으로 무시 가능한 수준이다. 따라서 사이클 슬립이 발생하지 않은 데이터를 이용하면 IF 선형조합의 평균을 통해 수신기의 시계 drift 추정이 가능하다. 수신기의 시계 drift 추정에 앞서 사이클 슬립의 스크리닝 과정이 수행된다. 수신기 시계의 영향 없이, 사이클 슬립을 효과적으로 모니터링 하기 위해서 2차 차분 관측값\(( ^i ∇^j Δ\tilde{Φ}_{IF}^i)\)이 사용된다. 여기서 \(^i ∇^j\)은 위성 \(i\)와 \(j\) 사이의 Single difference를 의미한다. Eq. (6)을 만족하는 IF 선형조합만을 사용하여 수신기의 시계 drift를 계산하며, 여기서 사이클 슬립 탐지의 기준으로 \(3σ_{i∇jΔ\tilde{Φ}_{IF}^i}\)가 사용된다.

\(| ^i ∇^j Δ\tilde{Φ}_{IF}^i |<3\cdot σ_{i∇jΔ\tilde{Φ}_{IF}^i}=3\cdot2\cdot\sqrt{(a_1 σ_{Φ_1} )^2+(a_2 σ_{Φ_2} )^2}\)                                               (6)

\(σ_{c\cdotΔ\widehat{dt}}^2={1\over m} σ_{Δ\tilde{Φ}_{IF}^i}^2={1\over m} (a_1^2 σ_{∆Φ_1}^2+a_2^2 σ_{∆Φ_2}^2)\)                                                                     (7)

이렇게 최종 계산된 수신기의 시계 drift의 오차 수준은 오차 전파 법칙을 적용하면 Eq. (7)과 같이 표현할 수 있다. 최종적으로, Eq. (8)과 같이 추정된 수신기의 시계 drift를 Eq. (3)에서 차분하여 보정된 반송파의 잔차를 의미하는 Carrier Phase Residual (CPR)을 계산할 수 있다.

\(CPR_1^i=Δ\tilde{Φ}_1^i (t)-cΔ\widehat{dt}=λ_1 ΔN_1^i-ΔI+ΔT+ε_{PRC_1^i}\)                                              
\(CPR_2^i=Δ\tilde{Φ}_2^i (t)-cΔ\widehat{dt}=λ_2 ΔN_2^i-γΔI+ΔT+ε_{PRC_2^i}\)                                       (8)

3. CYCLE-SLIP DETECTION AND COMPENTATION APPROACHES

앞서서 기준국 수신기의 정확한 위치 및 방송궤도력을 기반으로 반송파 관측값에 포함된 각종 오차를 보정하고, 시간 차분을 수행하여 CPR을 계산하였다. 이렇게 계산된 CPR을 이용하여 사이클 슬립을 검출하고, 보상하는 일련의 과정을 소개한다. 본 알고리즘은 대한 보다 자세한 설명은 Kim et al. (2018)을 통해 확인 가능하다.

3.1 Cycle-Slip Detection Using the Ionospheric Acceleration

CPR을 이용하여 사이클 슬립을 탐지하기 위하여 두가지 전리층 기반 선형조합을 적용한다. 첫번째는 L1, L2 신호의 차분을 통해 공통으로 포함되어 있는 기하 항을 소거하는 GF 선형조합이다.

\(ΔI^-=b_1^-\cdot CPR_1^i+b_2^-\cdot CPR_2^i=ΔI+{1\over γ-1} (λ_1 ΔN_1^i-λ_2 ΔN_2^i )+ε_{ΔI^-}\)                         (9)

여기서 \(b_1^-={1\over γ-1}\)이고, \(b_2^-=-{1\over γ-1}\)이다. Eq. (9)에서 확인되는 바와 같이 GF 선형조합에는 주파수에 종속적인 전리층과 모호정수 항만이 포함되므로 이중 주파수 기반 사이클 슬립 검출 방법으로 널리 사용된다. Kim et al. (2018)은 GF 선형조합에서 L1과 L2의 사이클 슬립이 음(negative)의 방향으로 더해진다는 점에서 Ionosphere Negative (IN) 선형조합으로 명명하였다. 일반적으로 IN 선형조합의 사이클 슬립 항은 전리층의 변화 속도와 관측 오차보다 크므로 시간에 따른 GF 값을 모니터링하여 사이클 슬립을 검출할 수 있다. 그러나 IN 선형조합은 기본적으로 각 주파수의 파장이 곱해진 L1, L2 신호의 사이클 슬립 차이\((λ_1 ΔN_1^i-λ_2 ΔN_2^i)\)를 모니터링 하므로 L1, L2 신호에 각각 (77n, 60n) 크기의 사이클 슬립이 발생한 경우 모니터링 값을 0으로 만들어 사이클 슬립의 탐지가 불가능하다. 이 밖에도, (4, 3), (5, 4), (9, 7) 등의 사이클 슬립 조합은 IN 값이 0.005 m/s\(^2\) 이하의 작은 값으로 계산되어 탐지가 불가능하다. 따라서 IN 선형조합의 단점을 보완할 수 있는 두번째 선형조합의 사용이 불가피하다.

일반적으로 사이클 슬립 탐지를 위해 IN 선형조합과 함께 MW 선형조합이 널리 사용된다. MW 선형조합은 반송파 관측값에 더불어 상대적으로 정확도가 떨어지는 의사거리 관측값을 사용하여 작은 크기의 사이클 슬립의 경우 탐지에 어려움이 따른다. 이에 Kim et al. (2018)은 높은 정확도의 반송파 관측값을 기반으로 L1, L2 주파수의 사이클 슬립의 양(positive)의 방향 합을 모니터링하는 Ionosphere Positive (IP) 선형조합을 처음 제안하였다. Eq. (10)에서 확인되는 바와 같이 IP 선형조합에는 대류층 지연오차 항 \(ΔT\)가 제거되지 않고 남아 있기는 하지만, 대류층 지연 오차의 변화량은 관측 간격이 짧은 기준국 관측값(1초 이하)을 이용하는 경우 무시할만한 수준이라고 할 수 있다.

\(ΔI^+=b_1^+\cdot CPR_1^i+b_2^+\cdot CPR_2^i\)                                                                                                                          
\(=-ΔI+1/2 (1+1/γ)ΔT+1/2 (λ_1 ΔN_1^i+1/γ λ_2 ΔN_2^i )+ε_{ΔI^+}\)                                 (10)

여기서 \(b_1^+={1\over2}\)이고, \(b_2^+={1\over2γ}\)이다. L1, L2 주파수의 사이클 슬립 항의 양(positive)의 방향 합을 모니터링 하는 IP 선형조합은 음의 방향으로의 합을 모니터링 하는 IN 선형조합과 같이 사용하였을 때 보다 효과적으로 사이클 슬립 탐지가 가능하다. Fig. 1은 사이클 슬립 조합에 따른 IN과 IP 선형조합 모니터링 값의 민감도 패턴을 보여준다. 사이클 슬립이 발생하지 않았음을 의미하는 (0, 0)을 기준으로 IN 선형조합의 민감도 패턴은 L1, L2에 같은 방향의 사이클 슬립이 발생하는 경우 모니터링 값이 작아지는 패턴을 보인 반면, IP 선형조합은 L1, L2에 반대 방향의 사이클 슬립이 발생하였을 때 모니터링 값이 작아지는 패턴을 보인다. 즉, IN, IP 선형조합은 서로 반대되는 민감도 패턴을 보이며, 따라서 두 선형조합을 상호 보완적으로 사용하면 사이클 슬립을 보다 효과적으로 탐지할 수 있다.

Fig. 1. IN and IP combinations’ sensitivity patterns depending on GPS L1 and L2 Cycle-Slip pairs.

앞서 전리층 폭풍이 발생하지 않은 상황에서 전리층 지연값의 변화 속도는 사이클 슬립에 의한 영향과 비교하여 무시할만한 수준이라고 언급한 바 있다. 그러나 Olynik et al. (2002)과 Olynik (2002)에 따르면 전리층 폭풍이 발생하는 경우 전리층 지연값의 변화율은 4배 이상 커질 수 있다. 따라서 전리층 폭풍이 발생하는 경우에도 사이클 슬립을 효과적으로 탐지하기 위해서 IN, IP 선형조합의 2차 시간 차분을 수행한다.

\(Δ^2 I^-=Δ^2 I+{1\over γ-1} (λ_1 ΔN_1^i-λ_2 ΔN_2^i )+ε_{Δ^2 I^-}\)                                                                                
\(Δ^2 I^+=-Δ^2 I+{1\over 2} \left(1+{1\over γ}\right) Δ^2 T+{1\over 2} \left(λ_1 ΔN_1^i+{1\over γ} λ_2 ΔN_2^i \right)+ε_{Δ^2 I^+}\)                         (11)

여기서 \(Δ^2 I^-\)와 \(Δ^2 I^+\)는 각각 IN과 IP 선형조합의 2차 시간 차분을 통해 얻어진 전리층 가속도이며, 전리층 가속도는 전리층 폭풍이 발생한 상황에서도 무시할만한 작은 값으로 나타난다. 또한, 시간에 따라 천천히 변화하는 대류층의 경우는 2차 시간 차분을 통해 더욱 그 영향이 작아지므로 사이클 슬립이 발생하지 않은 경우에는 IN과 IP의 전리층 가속도 값은 관측 오차만 포함되게 된다. 따라서, IN과 IP의 전리층 가속도 값을 사이클 슬립을 탐지하기 위한 모니터랑 값으로 정하였으며, Eq. (12)의 임계값을 초과하는 경우 사이클 슬립이 발생하였다고 할 수 있다.

\(|Δ^2 I^- |> T^-=3\cdot\max⁡ (σ_{Δ^2 I^-} )\)                                                                                  
\(|Δ^2 I^+ |> T^+=3\cdot\max⁡ (σ_{Δ^2 I^+} )\)                                                                           (12)

여기서 \(T^-\)와 \(T^+\)는 각 모니터링 값의 임계값이며, \(\max⁡(σ_{Δ^2 I^-})\)와 \(\max⁡( σ_{Δ^2 I^+})\)는 각각 IN과 IP 선형조합으로 얻어진 모니터링 값의 최대 오차 수준을 의미한다. 최대 오차 수준의 적용은 사이클 슬립이 누락되는 확률을 최소화하기 위한 것이다.

3.2 Error Propagation for the Threshold Determination

다음은 오차 전파 법칙을 기반으로 Eq. (12)에서 임계값 설정 시 사용한 모니터링 값의 오차 수준을 결정하는 일련의 과정을 설명한다. Eqs. (9-10)에서 IN, IP 선형조합을 적용하여 계산한 전리층 변화율은 Eq. (13)과 같이 일반화하여 표현 가능하다.

\(ΔI_{(b1,b2)}=b_1 (Δ\hat{Φ}_1-cΔ\widetilde{dt})+b_2 (Δ\hat{Φ}_2-cΔ\widetilde{dt})\)                                   (13)

일반화된 Eq. (13)에 오차 전파 법칙을 적용하고 Eq. (7)의 수신기의 시계 drift의 분산을 적용하면 전리층 변화율의 오차 수준을 Eq. (14)와 같이 표현할 수 있다.

\(σ_{ΔI_{(b1,b2)}}^2=b_1^2 σ_{Δ\hat{Φ}_1}^2+b_2^2 σ_{Δ\hat{Φ}_2}^2+{(b_1+b_2 )^2\over m} (a_1^2 σ_{ΔΦ_1}^2+a_2^2 σ_{ΔΦ_2}^2 )\)                                   (14)

여기서 \(σ_{Δ\hat{Φ}_1}^2\)와 \(σ_{Δ\hat{Φ}_2}^2\)는 각각 GPS L1, L2 반송파 관측값 변화율의 오차 수준을 의미한다. 또한, \(m\)은 수신기 시계 drift 계산시 사용 가능한 관측값의 수를 의미하는데, \(m\)이 1이 되었을 때 오차 수준이 최대가 된다. 따라서, IN, IP 선형조합 전리층 변화율의 최대 오차 수준은 Eq. (15)와 같이 표현 가능하다.

\(\max⁡(σ_{ΔI_{(b1,b2)}}^2 )=(b_1^2+a_1^2 (b_1+b_2 )^2 ) σ_{ΔΦ_1}^2+(b_2^2+a_2^2 (b_1+b_2 )^2 ) σ_{ΔΦ_2}^2\)                                   (15)

\(\max⁡(σ_{Δ^2 I_{(b1,b2)}}^2 )=6(b_1^2+a_1^2 (b_1+b_2 )^2 ) σ_{Φ_1}^2+6(b_2^2+a_2^2 (b_1+b_2 )^2 ) σ_{Φ_2}^2\)                                   (16)

이어서 2차 시간 차분 값의 분산은 6배 증가\((σ_{Δ^2 Φ}^2=6\cdot σ_Φ^2)\)된다는 사실을 고려하여 사이클 슬립 탐지를 위한 최종 모니터링 값인 전리층 가속도의 오차 수준은 Eq. (16)과 같이 계산된다. Eq. (16)에서 \(σ_{Φ_1}\)과 \(σ_{Φ_2}\)는 GPS L1, L2 신호의 오차 수준을 의미하며, 각각 3 mm와 3.85 mm \(\left({λ_2\over λ_1} \cdot σ_{Φ_1}=3.85~ {\rm mm}\right)\)를 적용할 경우 최종 임계값은 Eq. (17)과 같다.

\(T^-=3\cdot\max⁡ (σ_{Δ^2 I^-}^2)=0.055~ {\rm m/s}^2\)                                                                        
\(T^+=3\cdot\max⁡ (σ_{Δ^2 I^+}^2)=0.059~{\rm m/s}^2\)                                                                 (17)

3.3 Cycle-Slip Identification by Combining the IP and IN Combinations

앞서 설명한 바와 같이 서로 다른 특성을 갖는 IN, IP 선형조합을 이용하면 각 선형조합의 단점을 상호 보완하여 교차 검증이 가능하므로 효과적으로 사이클 슬립을 탐지할 수 있다. L1, L2 주파수 데이터에 포함된 정확한 사이클 슬립의 크기를 결정하기 위하여 다음의 과정이 수행된다. IN, IP 선형조합을 이용하여 독립적으로 탐지된 사이클 슬립은 Eq. (18)과 같이 표현 가능하다.

\(y=\left[\begin{matrix}Δ^2 I^-\\ Δ^2 I^+ \end{matrix}\right]=\left[\begin{matrix}b_1^- λ_1&b_2^- λ_2 \\ b_1^+ λ_1&b_2^+ λ_2 \end{matrix}\right]\left[\begin{matrix}∆N_1\\ ∆N_2 \end{matrix}\right]+\left[\begin{matrix}ε_{Δ^2 I^-}\\ ε_{Δ^2 I^+} \end{matrix}\right]=Ax+ε,\)                                  
\(ε\sim N (0,Q_{yy} ),~~~ Q_yy={\rm diag} (σ_{Δ^2 I^-}^2, σ_{Δ^2 I^+}^2 )\)                                                       (18)

여기에서 \(y\)는 앞서서 사이클 슬립으로 검출된 전리층 가속도 값을 포함하는 관측 벡터이며, \(x\)는 미지수 벡터로 L1, L2 주파수에 포함된 사이클 슬립을 포함한다. 또한, \(A\)는 관측 벡터와 미지수 벡터의 관계를 보여주는 디자인 행렬이며, \(Q_{yy}\)는 관측 벡터의 공분산 행렬이다. 각 주파수에 포함된 사이클 슬립의 크기는 다음과 같이 가중 최소 제곱법을 적용하여 추정 가능하다.

\(\hat{x}=\left[\begin{matrix}∆\widehat{N_1} \\(∆\widehat{N_2} \end{matrix}\right]=(A^T Q_{yy}^(-1) A)^{-1} A^T Q_{yy}^{-1} y\),                                                       (19)

Eq. (19)를 이용하여 추정된 각 주파수의 사이클 슬립은 부동소수점 값을 가지므로, 다음과 같이 부동소수점 값을 반올림하여 정수형 사이클 슬립으로 고정(fixed)한다.

\(\left[\begin{matrix}\widehat{{\rm CS}_1} \\ \widehat{{\rm CS}_2} \end{matrix}\right]=\left[\begin{matrix}{\rm round} (\widehat{∆N_1} )\\ {\rm round} (\widehat{∆N_1} ) \end{matrix}\right]\)                                                                           (20)

정수 사이클 슬립의 추정 후, 단순 이상 현상과 사이클 슬립의 영향을 분별하기 위한 검증 과정이 추가적으로 수행된다. 사이클 슬립 발생 외에 또다른 이상으로 인해 특정 epoch에 갑작스러운 오차가 포함되는 경우가 발생할 수 있는데, 이러한 경우 시간 차분 관측값에서 사이클 슬립과 유사한 점프 현상으로 나타날 수 있다. 검증을 위해서 반송파 관측값에 Eq. (20)에서 추정된 정수 사이클 슬립을 보상한다. 이어서, 사이클 슬립이 보상된 반송파 데이터를 사용하였을 때, IN, IP 모니터링 값이 임계값 내에 들어오는지 확인한다. IN와 IP 두가지 선형조합에서 모두 임계값 이내에 들어오는 경우에만 최종 사이클 슬립으로 판단하며, 그렇지 않은 경우에는 단순 이상값(outlier)으로 판별한다.

\(|Δ^2 I_{repaired}^- |<T^-=3σ_{Δ^2 I^-}=0.055~{\rm m/s}^2\)                                                                    
\(|Δ^2 I_{repaired}^+ |<T^+=3σ_{Δ^2 I^+}=0.059~{\rm m/s}^2\)                                                             (21)

이렇게 반송파 관측값을 기반으로 기준국의 무결성을 위한 사이클 슬립 탐지 및 판별을 수행하였다. 이러한 일련의 과정은 Fig. 2에서 확인 가능하다.

Fig. 2. Block diagram of the Cycle-Slip detection and identification algorithm for the integrity monitoring of the reference station.

4. Experimental validation

4.1 Simulation Environment

기준국 무결성 감시를 위해 본 논문에서 적용한 반송파 데이터 기반 사이클 슬립 탐지 알고리즘의 성능 평가를 위하여 제주(JEJU) 상시관측소에서 Trimble의 NetR8 수신기와 Zephyr Geodetic Model 2 안테나를 사용하여 1초 간격으로 수집된 데이터를 사용하였다. 특히, 본 실험에서는 GNSS의 각종 오차가 증가하는 특성으로 인해 사이클 슬립 검출이 까다롭다는 점을 고려하여 저앙각에서 수신된 GPS PRN 21, 22, 30 위성의 L1, L2 이중 주파수의 반송파 데이터를 사용하였다. 또한, 본 논문에서 적용한 사이클 슬립 알고리즘이 전리층 비활동기뿐만 아니라 전리층 지연의 값이 다이나믹하게 변화하는 전리층 활동기 동안에도 잘 작동하는지 확인하기 위하여 2017년 9월 8일과 9일 관측된 데이터를 이용하였다. 2017년 9월 8일은 강한 지자기폭풍이 발생한 날로, Fig. 3에서 확인되는 바와 같이 Kp 지수가 최대 8.333로 확인되었는데, 여기서 Kp 지수는 지구 자기장이 태양 입자와 만나 만드는 교란의 정도를 0~9등급으로 나타내는 지수이다. 반면 2017년 9월 9일은 Kp 지수가 3 이하로 확인되었으며, 통상 Kp 지수 3 이하는 지자기 활동이 활발하지 않은 안정적인 상황을 의미한다.

Fig. 3. The geomagnetic indices during 8-9 September, 2017.

사이클 슬립 검출 성능을 평가하기 위해서 저앙각에서 수신된 GPS PRN 21, 22, 30 위성의 L1, L2 반송파 관측값에 탐지가 까다로운 사이클 슬립 조합을 500 epoch 간격으로 임의로 추가하였다. 탐지가 까다로운 사이클 슬립 조합의 선정을 위해서 기타 전리층, 대류권, 관측 오차 등이 포함되지 않았음을 가정하고, GPS L1, L2 반송파에 각각 -10에서 10 cycle의 사이클 슬립 발생에 따른 IN, IP 선형조합의 이론적 모니터링 값을 Eq. (22)와 같이 계산하였다.

\(MV_{IN}={1\over γ-1} (λ_1 ΔN_1^i-λ_2 ΔN_2^i )\)                                                                          
\(MV_{IP}={1\over 2} \left(λ_1 ΔN_1^i+{1\over γ} λ_2 ΔN_2^i \right)\)                                                                 (22)

이어서 해당 모니터링 값이 각 임계값 이내로 계산되는 경우를 탐지가 까다로운 사이클 슬립 조합으로 선정하였다. 이렇게 최종 선정된 15개의 사이클 슬립 조합은 저앙각에서 수신된 GPS PRN 21, 22, 30 위성의 L1, L2 반송파 관측값에 500 epoch 간격으로 추가되었다. Table 1에 앞서 선정된 사이클 슬립 조합을 적용하였을 때, IN, IP 선형조합의 이론적 모니터링 값을 정리하였다. Table 1에서 정리된 바와 같이 (9, 7) 사이클 슬립 조합은 IN 선형조합의 모니터링 값이 0.49 cm/s\(^2\)의 작은 값으로 Eq. (17)의 IN 선형조합의 임계값으로는 탐지가 불가능할 것으로 예상된다. 이와는 반대로 (1, -2) 사이클 슬립 조합은 IP 선형조합의 모니터링 값이 -5.31 cm/s\(^2\)로 계산되어 Eq. (17)의 IP 선형조합의 임계값을 적용하였을 때, 탐지가 불가능할 것으로 예상된다.

Table 1. Insensitive Cycle-Slip pairs and their theoretical monitoring values.

Cycle-Slip pairs		Theoretical monitoring values
L1 [cycle]	L2 [cycle]	IN [cm]	IP [cm]
1 1 2 2 3 4 4 5 5 5 6 6 7 8 9	-2 -1 -3 -2 -4 -5 3 -7 -6 4 -8 -7 -9 -10 7	104.91 67.16 172.07 134.32 239.24 306.40 4.41 411.31 373.56 -3.92 478.47 440.72 545.63 612.80 0.49	-5.31 2.10 -3.21 4.20 -1.11 0.99 60.30 -4.32 3.09 77.23 -2.22 5.19 -0.12 1.98 137.53

 

4.2 Experimental Results

저앙각에서 관측된 GPS PRN 21, 22, 30 위성의 반송파 데이터에 매 500 epoch 마다 추가된 사이클 슬립 탐지 및 식별 결과를 분석하여 본 논문의 반송파 기반 사이클 슬립 검출 알고리즘의 성능을 평가하였다. 전리층 비 활동기 및 활동기 데이터를 이용한 사이클 슬립 탐지 결과는 각각 Figs. 4와 5의 위성별 IN, IP 모니터링 값의 시계열 그래프를 통해 확인 가능하다. 해당 그림에서 위쪽에 위치한 빨간색 시계열은 IN 선형조합 모니터링 값이며, 아래쪽의 파란색 그래프는 IP 선정조합의 모니터링 값이다. 추가로, 그래프에 회색 쉐도우 영역으로 표시된 임계값에 의해 탐지가 되지 않은 사이클 슬립을 검은색 '*'로 표시하였다. 전체적으로 분석 결과 그래프를 살펴보면, 전리층 활동 유무에 관계없이, 임의로 추가된 15개 사이클 슬립 조합 모두 IN, IP 선형조합 중 적어도 하나의 선형조합에 의해 분명히 탐지되었음을 확인할 수 있다. 따라서, IN, IP 두가지 선형조합을 상호 보완적으로 사용하는 알고리즘을 통해 GNSS의 각종 오차가 증가하는 저앙각에서 발생한 까다로운 사이클 슬립 조합도 효과적으로 탐지할 수 있었다. 이어서 L1, L2 주파수 관측값에 포함된 크기를 최종 결정하였다. 그 결과 전리층 활동 유무와 관계없이 임의로 추가한 15개 사이클 슬립 조합의 크기를 모두 성공적으로 결정할 수 있었다 (Table 2).

Fig. 4. CS detection results under the normal condition at JEJU station.

Fig. 5. CS detection results under the ionospheric storm at JEJU station.

Table 2. Cycle-Slip identification results under normal condition and ionospheric storm.

PRN	Inserted CS [cycle]	Ionospheric storm				Normal condition
		Ele [deg]	Estimated float CS [cycle]		Validated integer CS [cycle]	Ele [deg]	Estimated float CS [cycle]		Validated integer CS [cycle]
			\(∆N_1\)	\(∆N_2\)			\(∆N_1\)	\(∆N_2\)
G21	(1, -2) (1, -1) (2, -3) (2, -2) (3, -4)	7.15 9.10 11.02 12.89 14.69	1.02 1.02 1.99 2.02 3.04	-1.97 -0.99 -3.03 -1.98 -3.96	(1, -2) (1, -1) (2, -3) (2, -2) (3, -4)	8.46 10.39 12.28 14.10 15.84	1.00 1.02 2.00 2.04 3.00	-1.99 -1.01 -3.01 -1.97 -3.99	(1, -2) (1, -1) (2, -3) (2, -2) (3, -4)
G23	(4, -5) (4, 3) (5, -7) (5, -6) (5, 4)	7.33 9.48 11.67 13.87 16.06	4.00 4.00 5.00 4.99 5.01	-4.99 3.00 -6.99 -6.02 4.00	(4, -5) (4, 3) (5, -7) (5, -6) (5, 4)	7.27 9.43 11.61 13.81 16.00	3.99 4.00 4.97 5.01 4.99	-4.96 3.01 -7.02 -5.98 4.00	(4, -5) (4, 3) (5, -7) (5, -6) (5, 4)
G30	(6, -8) (6, -7) (7, -9) (8, -10) (9, 7)	7.23 9.33 11.44 13.55 15.64	6.05 6.00 6.95 8.01 9.02	-7.97 -7.00 -9.04 -9.99 7.02	(6, -8) (6, -7) (7, -9) (8, -10) (9, 7)	7.23 9.33 11.44 13.56 15.64	6.00 5.99 7.01 8.00 9.01	-8.00 -7.01 -9.00 -9.99 7.00	(6, -8) (6, -7) (7, -9) (8, -10) (9, 7)

 

5. CONCLUSIONS

사이클 슬립의 발생은 기준국 데이터의 무결성 유지에 심각한 영항을 미치므로 연속적이고 정확한 사이클 슬립의 감시 및 보상이 요구된다. 본 논문에서는 기준국 데이터의 무결성 감시에서 중요한 기준 값으로 사용되는 반송파 측정치의 무결성 감시를 위해 이중 주파수 반송파 데이터를 기반으로 하는 사이클 슬립 검출 및 식별 알고리즘을 제안하였다.

실시간, 고성능 사이클 슬립 검출 및 식별을 위해서 높은 정확도의 반송파 데이터만을 사용하는 IN, IP 선형조합을 적용하였으며, 이차 시간 차분된 값을 모니터링 값으로 사용하여 전리층 및 대류증 지연의 영향을 최소화하였다. 이를 통해 모니터링 값에 포함된 오차 수준을 낮추기 위한 추가 과정 없이도 오탐지를 최소화하였다. 또한, 사이클 슬립 조합에 따른 민감도 패턴이 서로 반대되게 나타나는 IN, IP 선형조합을 사용한 교차 검증을 적용하여 탐지가 불가능한 사이클 슬립 조합 없이 모든 사이클 슬립 검출을 가능하게 하였다. 사이클 슬립 탐지를 위한 임계값은 모니터링 값의 최대 오차 수준의 3배로 결정하였으며, 모니터링 값의 오차 수준은 오차 전파 공식을 적용하여 계산되었다. 최종적으로 GPS L1, L2 신호를 사용하였을 때의 임계값은, IN 선형조합의 경우는 0.055 m/s\(^2\)으로, IP 선형조합의 모니터링 값의 경우는 0.059 m/s\(^2\)으로 계산되었다. 각 주파수 데이터에 포함된 정확한 사이클 슬립의 크기를 결정하기 위해서는 가중 최소 제곱법을 적용하였다. 전리층 비활동기 및 활동기 데이터를 이용하여 논문에서 제시하는 기준국용 사이클 슬립 알고리즘의 성능 평가를 수행하였다. 각종 오차의 수준이 증가하는 저앙각 반송파 데이터에 검출이 까다로울 것으로 판단되는 15개 사이클 슬립 조합을 임의로 추가하고, 이렇게 추가된 사이클 슬립의 탐지 및 식별 여부를 확인하였다. 그 결과 전리층 활동 유무와 관계없이 임의로 추가된 모든 사이클 슬립 조합에 대하여 성공적으로 탐지 및 식별이 가능함을 확인할 수 있었다.

본 논문에서 제안하는 알고리즘이 실시간으로 고성능 감시를 수행해야 하는 기준국 데이터 무결성을 위한 이중 주파수 사이클 슬립 감시 및 식별에 적합한 알고리즘이라고 할 수 있다. 다만, 본 알고리즘을 기준국 무결성 감시 목적으로 사용하기 위해서는 장기간 관측된 데이터를 이용한 전리층 오차, 궤도/시계, 대류층 오차 등의 영향 분석을 진행하고, 이를 기반으로 하는 실질적인 임계값 설정이 필요할 것으로 생각된다.

AUTHOR CONTRIBUTIONS

Conceptualization, S.-K. Kim and D. Kim; methodology, S.-K. Kim and D. Kim; software, S.-K. Kim and S. C. Bu; validation, S.-K. Kim and B. Kim; formal analysis, S.-K. Kim; investigation, S.-K. Kim and D. Kim; resources, S.-K. Kim and S. C. Bu; data curation, S.-K. Kim and B. Kim; writing—original draft preparation, S.-K. Kim; writing—review and editing, S.-K. Kim, S. C. Bu and C. S. Lee; visualization, S.-K. Kim; Supervision, C. S. Lee.

CONFLICTS OF INTEREST

The authors declare no conflict of interest.