Accuracy Improvement of Multi-GNSS Kinematic PPP with EKF Smoother

Abstract

The extended Kalman filter (EKF) is widely used for global navigation satellite system (GNSS) applications. It is difficult to obtain precise positions with an EKF one-way (forward or backward) filter. In this paper, we propose an EKF smoother to improve the positioning accuracy by integrating forward and backward filters. For the EKF smoother experiment, we performed PPP using GNSS data received at the DAEJ reference station for a month. The effectiveness of the proposed approach is validated with multi-GNSS kinematic PPP experiments. The EKF smoother showed 35%, 6%, and 22% improvement in east, north, and up directions, respectively. In addition, accurate tropospheric zenith total delay (ZTD) values were calculated by a smoother. Therefore, the results from EKF smoother demonstrate that better accuracy of position can be achieved.

1. INTRODUCTION

전지구 위성항법시스템(Global Navigation Satellite System, GNSS)은 위치, 항법, 정밀시각 정보 제공 등에 활용되고 있다. 지상에서 정확한 3차원 위치정보를 얻기 위해서는 수신된 GNSS 관측자료를 정밀하게 처리하고 분석하는 작업이 필요하다. GNSS 관측자료를 이용한 위치결정 알고리즘은 최소자승법과 칼만필터가 대표적이다.

칼만필터(Kalman filter)는 1960년대 초에 개발된 알고리즘으로 현재 GNSS 항법, 제어 분야 등에서 널리 사용되고 있다. 칼만필터는 상태예측(State prediction)과 측정치 업데이트(Measurement update)를 반복적으로 수행하며 상태변수를 추정한다. 칼만필터는 시스템이 선형(Linear)이고, 모델이 가우시안(Gaussian) 분포를 따른다는 가정 하에서 사용할 수 있다. 그러나 실질적으로 많은 시스템이 비선형 구조를 가지고 있기 때문에 칼만필터를 적용하면 발산(Divergence) 등 여러 문제가 발생하는 경우가 있다. 따라서 칼만필터를 비선형에도 사용할 수 있도록 확장한 모델이 확장칼만필터(Extended Kalman filter, EKF)이다. EKF는 항법, 측위, 신호처리 등 GNSS 활용분야에 광범위하게 쓰이고 있다. Zhu et al. (2020)은 안드로이드 스마트폰으로 Multi-GNSS PPP를 수행했고, 자료처리 필터로 EKF를 사용했다고 보고했다. 또한 Erol (2020)은 동적환경에서 Multi-GNSS kinematic PPP의 성능을 비교하기 위해 EKF for-ward-backward-forward를 순차적으로 적용하였다. 그들은 초기 수렴되는 구간에서의 불안정성을 보완하기 위해 EKF forward와 backward 필터를 사용했다. 이처럼 EKF는 forward 및backward자료처리가 가능하기 때문에 forward의 단점들을 보완할 수가 있다 (Vaclavovic & Dousa 2015, Chiang et al. 2012).

Zhang et al. (2008)은 GPS와 dead reckoning를 이용한 차량의 통합항법을 위해 EKF smoother를 제안했고, EKF smoother가 EKF보다 더 높은 항법성능을 보인다는 결과를 제시했다. 또한 최근에 Zhang et al. (2020)은 GNSS PPP에 EKF forward와 back-ward 자료처리 결과를 통합한 항법결과를 제시했고, 모든 방향 성분별로 위치오차가 감소했다는 것을 발표했다.

GNSS PPP는 초기 모호수(Initial ambiguities)의 추정을 위해 수렴하는 시간이 필요하다. 따라서 EKF를 통해 추정변수들이 안정화 될때까지 최소 5분에서 30분까지의 시간이 필요하고, 때때로 1시간 이상의 수렴시간이 요구될 수도 있다. 초기 수렴시간을 줄이기 위한 방법으로는 Ionospheric-free 선형조합된 모호정수의 결정 (Ge et al. 2008, Geng et al. 2010), Multi-GNSS의 사용 (Cai & Gao 2013, Cai et al. 2015), 기상 수치 모델의 적용 (Wilgan & Geiger 2019), 그리고 고차항의 전리권 오차정보를 적용 (Hernandez-Pajares et al. 2007, Elmas et al. 2011) 등이 있다. 그러나 이 방법들은 초기 수렴되는 구간에서 추정변수의 불안정성을 완전하게 해결하지 못하고 있다.

본 연구에서는 Multi-GNSS kinematic PPP 자료처리시에 초기 추정변수의 불안정성을 보완하기 위해 EKF backward를 수행한다. 또한 EKF forward와 backward를 통합하는 EKF smoother 알고리즘을 구현한다. 그리고 EKF smoother와 EKF backward에 의해 계산된 위치오차를 서로 비교한다. 추가적으로 EKF smoother를 통해서 산출된 대류권 천정방향 지연오차(Tropospheric zenith total delay)를 국제 GNSS 서비스국(International GNSS Service, IGS)에서 제공하는 산출물과 비교 분석한다.

2. METHOD

2.1 EKF Smoother

EKF는 상태변수를 추정하는 최적화된 기법 중의 하나이다. EKF는 비선형 구조를 갖고 있는 시스템에 사용이 가능하도록 선형시스템에 사용하는 칼만필터를 확장한 개념이다. GNSS 측위 및 항법에서의 EKF는 표준화된 추정기법이라고 할 수 있다.

EKF 자료처리는 기본적으로 예측과 업데이트 부분으로 나뉜다. 예측 부분에서는 이전의 상태행렬과 제어행렬로부터 현시점의 상태변수들이 결정된다. 상태변수는 Choi & Yoon (2018)에 의해 정의된 것을 사용했다. 또한 예측된 상태변수들은 관측값을 적용하는 갱신 부분에 초기값으로 사용된다. 상태변수들을 갱신하기 위해서는 관측값, 칼만게인(Kalman Gain), 그리고 관측잡음이 고려된다 (Maatoug et al. 2017). 식 (1-4)는 EKF의 예측과 갱신을 나타낸 것이다.

예측 (Prediction):

\(\hat{x}_k=f(x ̃_{k-1}, \vec{u}_{k-1} )\)                                                                                                         (1)
\(\hat{P}_k=F_{k-1} \tilde{P}_{k-1} F_{k-1}^T+Q\)                                                                                                   (2)

갱신 (Update):

\(\tilde{x}_k= \hat{x}_k+K_k \hat{y}_k\)                                                                                                             (3)
\(\tilde{P}_k=(I-K_k H_k) \hat{P}_k\)                                                                                                         (4)

여기에서 \(f\)는 이전상태 \(x ̃_{k-1}\)에 대한 함수이고, \(\vec{u}_{k-1}\)는 제어벡터(control vector)이다. 만약 시스템이 외부 영향이 없다고 가정하면 제어벡터는 생략할 수 있다. \(P\)와 \(F\)는 각각 공분산 행렬(Covariance matrix)과 전이 행렬(Transition matrix)이고, \(Q\)는 시스템 잡음으로 가우시안 분포를 갖는다. \(K\)는 칼만게인, 그리고 \(y\)는 관측잔차이다. 또한 \(I\)와 \(H\)는 각각 정방행렬과 디자인 행렬을 의미한다.

위에서 설명한 것처럼 EKF에서 상태변수와 공분산 행렬은 시스템의 상태를 분석하는데 중요한 역할을 한다.

만약 어떠한 시스템이 작동을 \(k=1\)에 시작하여 특정 시점\((k=T)\)에 종료되었다고 가정하면, k가 1부터 \(T\) 시점까지 시간에 따라 진행을 하였기 때문에 ‘EKF forward’라고 한다. 반대로 종료시점(T)부터 초기시점 \(k=1\)까지 되돌아오면서 상태변수를 재추정하는 과정을 ‘EKF backward’라 한다. 본 연구에서 제안하는 EKF smoother는 EKF forward와 EKF backward에서 추정되는 상태변수와 공분산 행렬들이 결합되는 형태를 갖는다. 따라서 EKF smoother는 식 (5)와 (6) 처럼 표현된다 (Briers et al. 2010).

\(\tilde{P}_s= \left ( \tilde{P}_f^{-1}+\tilde{P}_b^{-1} \right )^{-1}\)                                                                                                         (5)

\(\tilde{x}_s=\tilde{P}_s \cdot \left ( \tilde{P}_f^{-1} \cdot \tilde{x}_f+\tilde{P}_b^{-1} \cdot \tilde{x}_b \right )\)                                                                                                 (6)

여기에서 \(\tilde{P}_f\)와 \(\tilde{P}_b\)는 각각 forward와 backward 과정에서 생성되는 공분산 행렬이다. 그리고 \(\tilde{x}_f\)와 \(\tilde{x}_b\)는 각각 forward와 backward과정의 상태변수를 의미한다. 최종적으로 EKF의 forward와 backward 과정에서 추정된 상태변수(\(\tilde{x}_f\)와 \(\tilde{x}_b\))와 공분산행렬(\(\tilde{P}_f\)와 \(\tilde{P}_b\)\(\))을 이용해서 평활화된 상태변수(smoothed state parameters), \(\tilde{x}_s\)가 계산된다.

2.2 GNSS PPP Strategy

본 연구에서는 GNSS PPP 자료처리에 EKF smoother를 적용하는 방법을 고안하였다. EKF smoother는 forward 자료처리 이외에 별도로 backward 과정이 필요하기 때문에 후처리 자료처리에 가깝다.

Table 1은 GNSS PPP 자료처리 방법들과 모델들을 기술한 것이다 (Choi & Lee 2021). GNSS PPP 자료처리에는 기준국에서 수신되는 모든 항법 신호가 고려된다. 여기에는 미국의 GPS L1/L2, 러시아의 GLONASS L1/L2, 유럽연합의 Galileo E1/E5a, 중국의 BeiDou B1/B3 그리고 일본의 QZSS L1/L5 이중주파수가 포함된다.

Table 1. Processing strategy for Multi-GNSS PPP (Choi & Lee 2021).

Item	Models / Methods
Estimator Observations Signals Elevation cutoff Sampling rate Satellite orbit and clock  Satellite/Receiver PCO & PCV Satellite differential code biases (DCB) Inter-system biases Phase wind-up Solid tide, ocean tide, pole tide Tropospheric error Mapping function Ambiguity	EKF forward/backward/smoother Un-differenced ionosphere free linear combination GPS L1/L2, GLONASS L1/L2, Galileo E1/E5a, BeiDou B1/B3, QZSS L1/L5 7o 300 sec GFZ products igs14.atx CODE monthly products Estimated as white noise, GPS as reference Wu et al. (1993) & nominal yaw-steering mode IERS conventions 2010 Zenith wet delay estimation with gradient parameters GMF/GPT2 Float

 

PPP에서 핵심이 되는 위성의 궤도와 위성시각 오차는 독일 지구물리연구소(GeoForschungsZentrum Potsdam, GFZ)에서 제공하는 산출물을 이용한다. 그리고 안테나에 대한 위상중심옵셋(phase center offset, PCO)과 위성중심변화(phase center variation, PCV)가 고려된다. 또한 C1 또는 C/A 코드를 P1코드로 변환하기 위해 항법위성의 하드웨어 바이어스(P1C1.DCB와 P1P2.DCB)를 스위스 베른대학의 Center for Orbit Determination in Europe (CODE)로부터 다운로드 한다.

항법위성과 지상의 수신기간의 상대적인 움직임에 의해 생성되는 위상말림효과(phase wind-up)는 Wu et al. (1993)이 고안한 방법을 적용한다. 그리고 지각, 해양, 극 조석에 의한 영향도 함께 고려한다.

GNSS PPP에서 추정되는 변수들 중의 하나는 천정방향의 대류권 지연오차(Zenith total delay)이다. 천정방향과 시선방향 오차 지연값의 전환을 위한 사상함수로는 Global Mapping Function (GMF)를 적용하고, 대류권 오차추정에 필요한 평균온도와 대기압력 값은 Global Pressure and Temperature 2 (GPT2) 모델에서 생성된 값을 이용한다. 그리고 Ionosphere free 선형조합된 반송파의 모호수(ambiguities)는 실수값(float values)으로 적용된다.

3. RESULTS AND ANALYSIS

DAEJ 국제 GNSS 기준국에서 2020년 1월 1일부터 1월 31일까지 수신된 GNSS 관측자료를 이용하여 Kinematic PPP 자료처리를 수행하였다. PPP 자료처리 과정에서 forward, backward 그리고 smoother 필터에 의한 측위결과를 서로 비교 하기위해 먼저 2020년 1월 1일에 수신된 자료를 처리하였다.

Fig. 1은 2020년 1월 1일 DAEJ 기준국에서 관측된 가시위성 수의 변화를 나타내고 있다. 수신된 항법위성은 GPS, GLONASS (GLO), Galileo (GAL), BeiDou (BDS) 그리고 QZSS (QZS)로 구성된다. GPS와 BDS 항법위성은 일부구간에서 최대 11기까지 그리고 GLO와 GAL은 최소 4기에서 최대 8기까지 수신이 되었다. 또한 QZS는 최소 2에서 최대 3기까지 수신이 되는 것을 알 수 있다. Fig. 1에서 볼 수 있듯이 검은색 실선은 전체 가시위성 수의 변화이고, 최소 26기에서 최대 36기까지 수신되었다. 물론 가시위성 수는 관측환경과 위성의 임계고도각(cutoff angle) 설정에 따라 달라질 수 있다.

Fig. 1. Changes of the number of visible satellites at DAEJ GNSS reference station on January 1, 2020.

Fig. 2는 DAEJ 기준국에서 모든 항법위성군의 조합에 의해 계산된 Dilution Of Precision (DOP)의 변화를 나타낸 것이다. DOP는 항법위성군의 기하학적인 배치 정도를 나타내는 지표로써 본 연구에서는 geometry DOP (GDOP)와 position DOP (PDOP)의 변화를 제시했다. DOP는 값이 작을수록 정밀도의 확률이 높아지는 것을 의미한다. Fig. 2에서 볼 수 있듯이 GDOP와 PDOP 값들은 모두 1.0 이하이고, 0.5와 1.0 사이에서 안정적인 것을 알 수 있다. Li et al. (2015)에 따르면 4개의 서로 다른 항법시스템의 조합으로 가시위성 수가 22에서 30기까지 바뀔 때, PDOP 값은 1.5 이하의 값을 갖는다고 보고한 바 있다. 우리의 결과는 관측된 위성 수가 그들의 결과에 비해 상대적으로 많고, QZS 항법시스템이 추가되었기 때문에 PDOP가 작게 계산된 것으로 판단된다.

Fig. 2. DOP changes at DAEJ station with multi-GNSS on January 1, 2020. The red and blue solid lines represent the GDOP and PDOP, respectively.

Fig. 3은 2020년 1월 1일 DAEJ 기준국에서 관측된 자료를 이용한 PPP 자료처리 결과이다. EKF의 forward, backward 그리고 smoother에 의해서 정밀 위치해가 각각 계산되었다. 본 연구에서는 정밀 위치해가 매 5분 간격으로 결정되어 저장되도록 설정했다. Fig. 3에서 회색실선은 forward 방법에 의해 결정된 위치오차의 시계열 정보이다. 여기에서 정밀위치해의 검증을 위한 참값은 국제 GNSS 서비스 기구에서 Software INdependent EXchange (SINEX) 형식으로 제공하는 산출물을 이용하였다. forward 방법은 초기 모호수들(Initial ambiguities)이 수렴하는 과정이 필요하기 때문에 위치해의 정확도가 떨어지는 단점이 있다. 반면 검은색 점선으로 표기된 backward 방법은 자료처리 시간을 되돌려서 위치해를 계산하기 때문에 초기 위치해의 정확도가 떨어지는 것을 보완하는 역할을 한다. smoother에 의한 정밀 위치해는 붉은색 실선으로 나타내었다. Fig. 3에서 알 수 있듯이 smoother가 forward와 backward에 비해 상대적으로 안정적인 결과를 산출한다는 것을 알 수 있다.

Fig. 3. GNSS PPP errors with three different methods (forward, backward and smoother) at DAEJ station on January 1, 2020. The solid gray line and black dashed lines represent the posi-tion errors with forward and backward, respectively. In addition, the results for the smoother filter are shown as the solid red lines. The gray solid and black dashed arrows indicate the direc-tion in GNSS data processing.

Table 2는 Fig. 3에서 제시된 위치오차를 각 성분별 root mean squares (RMS)값을 제시한 것이다. RMS값은 95%의 신뢰구간(confidence level) 내에서 계산되었다. forward 방법은 동서(E), 남북(N), 고도(U) 방향의 RMS값이 각각 2.28, 1.27, 그리고 2.06 cm로 산출되었다. 그리고 backward 방법은 RMS값이 각각 0.68, 0.47, 그리고 1.38 cm를 보였다. 결과적으로 backward가 forward에 비해 성능이 우수한 것으로 나타났다. 특히 동서와 남북방향 성분의 RMS값이 현격히 줄어드는 결과를 보이고 있다. smoother는 각 성분의 RMS값이 0.59, 0.45 그리고 1.10 cm로 산출되었다. 이러한 결과는 forward와 backward에 비해 위치해의 성능이 우수한 것을 알 수 있다. 특히 smoother는 backward에 비해 동서와 고도방향에서 좀 더 성능이 향상되는 특징을 보인다.

Table 2. Position RMS errors at DAEJ station with EKF forward, backward, and smoother on January 1, 2020 (unit: cm).

Method	RMS
	E	N	U
Forward Backward Smoother	2.28 0.68 0.59	1.27 0.46 0.45	2.06 1.38 1.10

 

앞서 우리는 1일 관측자료를 이용하여 세가지 방법을 통한 PPP의 위치해를 분석한 바 있다. smoother 방법에 대해 좀 더 신뢰성을 확보하기 위해서 2020년 1월 1일부터 31일까지(총 31일간)의 자료처리를 수행하였다.

Fig. 4는 2020년 1월 1일부터 31일까지 backward와 smoother에 의한 위치오차를 시계열로 나타낸 것이다. backward와 smoother에 의한 위치해는 각각 검은색 점선과 붉은색 실선으로 표기했다. Fig. 4에서 두 개의 방법을 서로 비교했을 때, smoother가 backward에 비해 모든 위치성분에서 안정적인 결과를 산출한다는 것을 알 수 있었다.

Fig. 4. PPP position errors at DAEJ station with backward and smoother methods from DOY (day of year) 1 to 31, 2020. The dotted line and solid red lines represent the backward and smoother solutions, respectively.

Table 3은 backward와 smoother에 의한 위치오차를 각 성분별 RMS값으로 나타낸 것이다. backward는 각 성분별 RMS값이 각각 0.72, 0.58, 그리고 1.47 cm로 계산되었다. 반면 smoother는 0.47, 0.54, 1.14 cm로 계산되어 backward에 비해 동서, 남북, 그리고 고도 방향에서 각각 35%, 6%, 22%의 향상된 결과를 보였다. 또한 smoother는 동서(E)와 고도(U)방향에서 위치정확도의 향상이 더 크게 나타났고, 남북방향 성분은 차이가 크지 않았다.

Table 3. Position RMS errors at DAEJ station with EKF backward and smoother from DOY 1 to 31, 2020 (unit: cm).

Method	RMS
	E	N	U
Backward Smoother	0.72 0.49	0.58 0.54	1.47 1.14

 

GNSS PPP에서 대류권 지연오차는 위치오차에 큰 영향을 준다. 특히 고도방향 성분의 위치오차와 직접적으로 관련이 있다 (Haase et al. 2003). PPP 자료처리를 통해서 대류권의 천정방향 지연오차(zenith total delay, ZTD)를 정밀하게 추정할 수 있다. Fig. 5는 2020년 1월 1일 하루동안 DAEJ 기준국에서의 대류권 ZTD의 변화를 보여주고 있다. forward와 backward에 의한 대류권 ZTD는 각각 회색 실선과 검은색 점선으로 표시했다. Fig. 5에서 forward에 의해 산출된 ZTD는 초기 모호수가 수렴되는 일부구간에서 불안정하게 산출되는 것을 알 수 있다. 그리고 화살표는 자료처리의 방향성을 의미한다. 붉은색 실선은 smoother에 의해 산출된 대류권 ZTD를 나타내고 있다. smoother에 의해 추정된 대류권 ZTD의 검증을 위해 IGS에서 제공하는 ZTD 산출물과 비교하였다. IGS는 하루단위로 ZTD 산출물을 제공하고 있고, ZTD는 매 5분 간격으로 생성된다.

Fig. 5. Tropospheric ZTD at DAEJ station with three different methods (forward, backward, and smoother) on January 1, 2020. The solid gray line and black dashed lines represent the ZTD changes with forward and backward, respectively. In addition, the results for the smoother filter are shown as the solid red lines. The gray solid and black dashed arrows indicate the direc-tion in GNSS data processing.

Fig. 6은 대류권의 천정방향 습윤 지연오차(zenith wet delay, ZWD)에 대한 공분산 값의 변화를 나타낸 것이다. GNSS PPP에서는 대류권 ZWD를 정밀추정하는 것이 매우 중요하다. 따라서 각 자료처리 필터에 따른 ZWD의 공분산 변화를 분석할 수 있다. Fig. 6에서 회색과 검은색 실선은 각각 forward와 backward에 의해 계산된 대류권 ZWD의 공분산을 의미한다. forward는 초기에 공분산 값이 안정화되는 시간이 필요하다. 또한 붉은색 실선은 smoother에 의해 산출된 공분산 값으로 앞서 forward와 backward에 비해 공분산 값이 상대적으로 작은 것을 알 수 있다.

Fig. 6. Covariances for tropospheric ZWD at DAEJ station with three different methods (for-ward, backward, and smoother) on January 1, 2020. The solid gray line and black lines repre-sent the covriances values with forward and backward, respectively. In addition, the results for the smoother filter are shown as the solid red line.

Fig. 7은 2020년 1월 1일부터 1월 31일까지 DAEJ 기준국에서의 대류권 ZTD 시계열을 나타낸 것이다. 검은색 점선은 IGS ZTD 값을 의미하고, 붉은색 점들과 실선으로 표기된 것은 본 연구에서 제안한 smoother에 의해 추정된 ZTD값을 의미한다. IGS ZTD와 smoother ZTD값의 변화는 매우 유사한 것을 알 수 있다. IGS와 smoother간에 RMS값과 표준편차는 각각 약 5.0 mm와 4.7 mm로 계산되었다. ZTD는 일일단위로 생성되기 때문에 일부 날짜 경계에서 불연속적인 특징을 보이기도 한다. 또한 IGS ZTD는 일부 순간적인 큰 변화가 나타나기도 한다. 그럼에도 불구하고 두 비교간의 RMS값이 5.0 mm로 산출된 것은 smoother에 의해 추정된 ZTD가 안정적이고 신뢰할 만하다는 것을 의미한다.

Fig. 7. ZTD time series at the DAEJ station from DOY 1 to 31, 2020. The dotted lines represent the IGS ZTD values. The red dotted and solid lines denote the smoother ZTD values.

4. CONCLUSIONS

본 연구에서는 GNSS PPP 자료처리에 EKF smoother를 적용하여 정밀 위치해와 대류권 ZTD를 산출하였다. 먼저 2020년 1월 1일부터 31일까지 DAEJ 기준국의 관측자료를 이용해서 Multi-GNSS kinematic PPP를 수행했다. smoother에 의해 산출된 위치해의 검증을 위해 IGS SINEX 참값과 비교하였고, 동서, 남북 그리고 고도방향의 RMS값이 각각 0.47, 0.54, 1.14 cm로 산출되었다. 또한 smoother는 backward와 비교하여 동서와 고도방향 성분에서 각각 35%와 22%의 위치정확도가 향상되었다.

추가적으로 GNSS PPP에서 smoother의 추정변수로써 대류권 ZTD를 산출하였고, IGS에서 제공하는 ZTD 산출물과 비교를 하였다. 본 연구에서 제안한 smoother ZTD와 IGS ZTD간의 RMS와 표준편차는 각각 5.0 mm와 4.7 mm 수준으로 신뢰할 만한 결과를 얻었다.

본 연구를 통해 우리는 EKF smoother가 Multi-GNSS kinematic PPP의 위치정확도 향상과 신뢰할 만한 대류권 ZTD를 산출한다는 것을 보였다. 따라서 EKF smoother는 준 실시간 또는 후처리 자료처리가 필요한 많은 비선형시스템에 적용이 가능할 것으로 판단된다.