Evaluation of RTK Methods for Moving Vehicles and Practical Recommendations

Abstract

Global Navigation Satellite Systems (GNSS) based precise positioning using Real Time Kinematic (RTK) technique has been proposed as an enabler of the formation operation of moving vehicles. In RTK methods, the integer ambiguity of GNSS carrier phase measurements must be resolved. Although there have been many proposed algorithms for the integer ambiguity resolution, the widelane combination of carrier phase measurements and LAMBDA methods have gained the most popularity in literatures when dual frequency GNSS measurements were used. In this paper, we evaluated five alternative methods to determine relative positions of moving base and rover receivers; the round-off scheme of widelane carrier phase, instant least-squares and Kalman filter-based LAMBDA with widelane carrier phase, instant least-squares and Kalman filter-based LAMBDA with dual frequency measurements. The paper presented the performance of each method using flight test data, which showed their strength and weakness in the aspects of time-to-first-fix, ambiguity resolution success ratio, and relative position errors. Based on that, we provided practical recommendations of RTK operations for moving vehicles.

1. 서론

전지구 위성 항법 시스템인 Global Navigation Satellite Systems (GNSS)는 정확하고 지속적인 위치, 항법 및 타이밍 (PNT) 서비스를 제공하기 위해 광범위하게 사용된다. 측량, 건설기계 제어, 정밀농업, 자율 주행, 건물 변위 측위 등을 포함한 정밀한 위치 정확도가 요구되는 사례에 위성항법 시스템은 많이 사용되고 있으며 (Olsen et al. 1999, Peyret et al. 2000), 센티미터 수준의 정확도를 제공하기 위해 Real-Time Kinematics (RTK) 또는 Precise Point Positioning (PPP) 기법이 사용된다. RTK는 이동국의 의사 거리 측정치만 사용하는 단독 측위와 달리 이동국과 기준국의 위성항법 신호 반송파 위상 측정치를 이용하여 센티미터 수준의 정밀한 상대 벡터를 계산한다. 반송파 위상 측정치를 위치결정에 사용하기 위해서는 우선 반송파 위상에 포함되어 있는 미지정수해를 결정해야 한다 (Montenbruck et al. 2002, Odolinski & Teunissen 2016). 이와 관련하여 검색 범위 축소와 같은 정수 모호성 해결을 돕기 위한 연구가 많이 행해졌으며, 미지정수 결정에 대표적인 기법으로는 Least-squares AMBiguity Decorrelation Adjustment (LAMBDA)가 있다 (De Jonge & Tiberius 1996, Teunissen 2006, Li et al. 2013a). 한편 PPP는 측정값 이외에 정밀한 위성 궤도력 정보, 대류 및 전리층 지연, 위성시계 보정정보 그리고 지각변동 보정모델등을 적용하여 높은 수준의 정확도를 제공하지만 대체로 위치측정 수렴 속도가 수분에서 수십분에 달할수 있으며 PPP 서비스가 제공되는 지역 이외에서는 사용할 수 없다는 한계를 가진다 (Zumberge et al. 1997, Bisnath & Gao 2009, Li et al. 2013b).

RTK는 이동국과 기준국의 위성항법 신호 코드 및 반송파 위상 측정치를 이중 차분하여 전리층 및 대류층 지연으로 인한 오차와 인공위성의 시계오차를 제거하므로 높은 정확도의 상대벡터 값을 제공할 수 있다. 이때, 상대벡터의 정확도는 반송파 위상 이중 차분 미지정수 추정 성능에 따라 달라지며, 이동국과 기준국 간의 기저선 거리가 멀어질수록 전리층 및 대류층에 대한 오차들의 상관관계가 적어지기 때문에 10∼20 km 이내의 기저선 (short baseline) 이 존재할 경우에 신뢰성 높은 미지정수해와 상대벡터를 계산할 수 있다 (Zhao et al. 2014, He et al. 2014, Odolinski et al. 2015).

일반적인 RTK 기법은 이동국과 고정된 기준국이 측정값을 동시에 받고 이동국의 위치를 계산할 때 기준국의 측정값으로부터 시각 보정된 위치를 계산한다. 기준국이 측정값을 송신하는 시간과 이동국이 측정값을 받는 시간의 기준국의 위치는 동일하므로 고정된 기준국의 위치는 이동국과 기준국의 상대벡터 계산에 영향을 미치지 않는다. 그러나 기준국과 이동국 사이의 거리가 멀수록 이동국의 위치오차가 증가하기 때문에 이동국의 정밀한 위치 측정을 위해 기준국이 이동하는 moving baseline RTK를 사용해야한다. 무인항공기의 편대 비행의 경우 각 기체의 정밀한 위치가 요구되는데, 이 경우 고정된 기준국 보다 임의의 기체 하나를 기준국으로 설정한 moving baseline RTK가 더 적합하다. 하지만 기존의 고정된 기준국을 사용했을 때와 비교해 비교적 정확한 기준국의 위치를 알 수 없기 때문에 정밀한 미지정수 및 상대벡터 추정에 비교적 어려움이 존재한다. 따라서 미지정수와 상대벡터를 계산하기 위한 상태변수를 기준국의 위치에 의존하지 않도록 구성하여, 이동국과 기준국의 시각 동기화 이후 기존 고정된 기준국 수신기를 이용하는 RTK 방식과 대체로 유사하게 구현할 수 있다. 올바른 미지정수해 도출을 위해 기준국 수신기의 위치가 정확히 측정된 후 RTK를 수행하는 방법도 제안되었으나 이는 초기 미지정수 결정에 긴 시간이 요구된다 (Psiaki & Mohiuddin 2005, Travis et al. 2008, D’Amico & Montenbruck 2010). Melbourne-Wübbena가 제시한 Melbourne-Wübbena 조합을 사용한 RTK 기법 또한 제안되었으며, 실제 Sea-based Joint Precision Approach and Landing 에도 사용된 바 있다 (Dogra et al. 2005). 하지만 LAMBDA 기법을 사용하여 미지정수를 결정하는 경우, 확장파장 조합을 측정치로 사용하는 것은 단순히 이중 주파수 측정치를 사용한 경우와 비교하여 이중 주파수 측정치의 일부 조합만을 탐색하는 것이기 때문에 (Teunissen 1997), 이중 주파수 측정치를 사용한 경우에 비해 확장파장 조합을 사용한 경우가 미지정수해 도출 성능이 우수하다 하기 어려우며, 상대벡터 도출시 증폭된 노이즈로 인해 계산된 상대벡터 값이 일반적인 RTK보다 정밀하지 않을 수 있다.

본 논문에서는 moving baseline RTK에 널리 사용되는 확장파장 기법과 고정된 기준국 위치에 의존하지 않는 최소제곱법 (least-squares) 및 칼만필터 기반 RTK 성능을 비교하고 저가의 수신기를 탑재한 무인항공기 편대 비행을 위해 고려될 사항을 제시하였다. 확장파장 조합은 반올림 (round-off) 기법과 최소제곱법 기반 LAMBDA기법, 칼만필터 기반 LAMBDA 기법을 적용하여 미지정수해를 도출하였고, 이중 주파수의 경우 확장파장 조합의 동일한 LAMBDA 기법을 사용하였다. 또한 기준국과 이동국 측정값의 시간차가 존재할 경우 동일한 GPS 시각을 기준으로 시각 동기화 (Time Synchronization)를 진행하였다. 미지정수의 정수해 검정 기법인 ratio test를 통과한 경우에만 LAMBDA의 고정 정수해 (fixed integer)를 사용하였고, 초기 위치 결정 시간 (Time-To-First-Fix, TTFF), 미지정수해 도출 성공 확률 및 기준국과 이동국 수신기 상대위치 벡터 오차의 제곱평균제곱근 (root mean square, RMS)을 성능지표로 사용하였다. 본 논문의 2장에서는 이중 주파수 위성 항법 반송파 위상 측정치 조합을 소개하고 이를 이용한 상대위치 벡터 산출 수식들을 나타내었으며, 3장에서는 moving baseline RTK 알고리즘검증에 사용된 데이터를 획득하기 위한 실험 환경 및 구성을 소개한 뒤, 4장에는 각 기법에 대한 상대항법 결과를 비교하였다.

2. 이중 주파수 위성 항법 반송파 위상 측정치를 이용한 상대 항법

이번 장에서는 이중 주파수 측정치 모델, 확장파장 조합 측정치 및 미지정수해 도출과정을 기술하였다.

2.1 이중 주파수 측정치 모델

위성 항법 L1과 L2 주파수 측정치들은 식 (1-4)와 같이 모델링 될 수 있다.

\(\Phi_{L1}=r-I_{L1}+T+(c\delta t_u-c\delta t^s )+\lambda_{L1} N_{L1}+\varepsilon_{\Phi_{L1}}\)                                                    (1)

\(\Phi_{L2}=r-I_{L2}+T+(c\delta t_u-c\delta t^s )+\lambda_{L2} N_{L2}+\varepsilon_{\Phi_{L2}}\)                                                  (2)

\(\rho_{L1}=r+I_{L1}+T+(c\delta t_u-c\delta t^s )+\varepsilon_{\rho_{L1}}\)                                                                          (3)

\(\rho_{L2}=r+I_{L2}+T+(c\delta t_u-c\delta t^s )+\varepsilon_{\rho_{L2}}\)                                                                          (4)

여기에서 \(r\)은 수신기로부터 위성까지의 실제 거리, \(I_{L1},I_{L2}\)는 L1과 L2의 전리층 지연 오차, \(T\)는 대류층 지연오차, \(c\delta t_u\)는 수신기 시계 오차, \(c\delta t^s\)는 위성 시계 오차, \( \lambda_{L1}, \lambda_{L2}\)는 L1과 L2의 파장, \(N_{L1},N_{L2}\)는 각각 L1과 L2 주파수의 반송파 위상 미지정수해이며 \(\varepsilon\)은 다중경로와 수신기 잡음에 의한 오차이다.

기준국과 이동국의 수신기 간의 단일 차분 (single difference, sd)에서는 기준국과 이동국의 거리가 10 km 이내인 경우 위성시계 오차 및 전리층과 대류층 지연 오차가 상쇄되며, 기준 위성과의 이중 차분 (double difference, dd)에서는 수신기 시계오차가 상쇄된다. 단기선 (short baseline)에서의 반송파 위상 이중 차분 측정치 \(\Phi_{dd}\)와 코드 이중 차분 측정치 \(\rho_{dd}\)는 식 (5-8)과 같다.

\(\Phi_{dd,L1}=-(e^i-e^j)x+\lambda_{L1} N_{dd,L1}+\varepsilon_{\Phi_{dd,L1}}\)                                                       (5)

\(\Phi_{dd,L2}=-(e^i-e^j)x+\lambda_{L2} N_{dd,L2}+\varepsilon_{\Phi_{dd,L2}}\)                                                       (6)

\(\rho_{dd,L1}=-(e^i-e^j)x+\varepsilon_{ρ_{dd,L1}}\)                                                                                   (7)

\(\rho_{dd,L2}=-(e^i-e^j)x+\varepsilon_{ρ_{dd,L2}}\)                                                                                    (8)

여기서 \(e^i\)는 i번째 위성과 이동국의 시선각 벡터, \(x\)는 기준국으로부터 이동국의 상대벡터, \(N_{dd,L1},N_{dd,L2}\)는 각각 L1과 L2의 이중 차분 미지정수이다.

미지정수 추정 시의 모호성은 반송파의 파장 길이에 따라 달라지며 파장이 길수록 추정치가 더 좋다 할 수 있다. 확장파장 조합의 파장은 L1과 L2 주파수의 파장에 비해 더 길기 때문에, 미지정수 추정 오차가 상대적으로 줄어든다. L1 의사거리를 사용하여 L1 반송파 위상의 미지정수를 추정할 때의 추정 오차 모호성은 \( \sigma(\hat{N}_{L1}) \approx 5 \) cycle인 반면 확장파장 조합의 미지정수 \(N_w\)를 narrowlane 의사거리로 측정할 때의 추정 오차는 \(\sigma(\hat{N}_w )\approx 1.2\) cycle로 추정 정확도가 증가하게 된다. 차분 전 반송파 위상 확장파장 조합 (Undifferenced phase widelane) \(\Phi_w\)과 코드 축소파장 조합 (Undifferenced code narrowlane) \(\rho_n\)은 식 (9), (10)과 같다.

\(\Phi_w={\lambda_L2 \over \lambda_{L2}-\lambda_{L1}} \Phi_{L1}-{\lambda_{L1} \over \lambda_{L2}-\lambda_{L1}} \Phi_{L2}=r-{\lambda_{L2}\over \lambda_{L1}}  I_{L1}+T+(c\delta t_u-c\delta t^s )+\lambda_w N_w+\varepsilon_{\Phi_w}\)                  (9)

\(\rho_n={\lambda_{L2} \over \lambda_{L2}\lambda_{L1}} \rho_{L1}+{\lambda_{L1} \over \lambda_{L2}+\lambda_{L1}} \rho_{L2}=r+{\lambda_{L2} \over \lambda_{L1}}  I_{L1}+T+(c\delta t_u-c\delta t^s )+\varepsilon_{\rho_n}\)                                         (10)

여기서 확장파장 조합의 파장은 \(\lambda_w={c \over f_{L1}-f_{L2}}={\lambda_1 \lambda_2 \over \lambda_2-\lambda_1}\)이며, \(N_w\)은 반송파 위상 확장파장 조합의 미지정수이다.

L1과 L2의 반송파 위상 확장파장 조합과 코드 축소파장 조합도 마찬가지로 단기선내에서 전리층과 대류층의 오차, 위성 시계와 수신기 시계 오차는 이중 차분을 통해 결과 상쇄된다. 이중 차분 반송파 위상 확장파장 조합 \(\Phi_{dd_w}\)과 이중 차분 코드 축소파장 조합 \(\rho_{dd_n}\)은 식 (11), (12)와 같다.

\(\Phi_{dd,w}=-(e^i-e^j)x+\lambda_w N_{dd,w}+ \varepsilon_{\Phi_{dd,w}}\)                                                             (11)

\(\rho_{dd,n}=-(e^i-e^j)x+\varepsilon_{\rho_{dd,n}}\)                                                                                     (12)

이중 차분 코드 축소파장 조합 공분산 \(R_{\rho_{dd,n}}\)과 이중 차분 반송파 위상 확장파장 조합 공분산 \(R_{\Phi_{dd,w}}\)은 각각 식 (13)과 (14)이며, 식 (15-24)의 과정으로 계산된다. 먼저 각 조합 측정치의 이동국과 기준국에 대한 단일 차분 측정 오차의 분산을 구하고, 단일-이중 차분 변환 행렬인 D 행렬을 곱하여 이중 차분된 측정 값들의 공분산을 계산한다.

\(R_{\rho_{dd,n}}=DR_{\rho_{sd,n}} D^T\)                                                                                         (13)

\(R_{\Phi_{dd,w}} = DR_{\Phi_{sd,w}} D^T\)                                                                                       (14)

여기서 \(R_{\rho_{sd,n}}\)와 \(R_{\Phi_{sd,w}}\)은 각각 축소파장 조합과 반송파 위상 확장파장 조합에서의 단일 차분 측정 오차 분산 행렬이고, 행렬 \(D\)는 GPS, GAL와 BDS 각 위성군에서 pivot 위성을 기준으로 단일 차분 형태에서 이중 차분 형태로의 변환 행렬이다. 만약 GPS 위성군에서 n개의 위성이 가시영역안에 있다면, 식 (16)과 같이 \(D^{GPS}\)는 \((n-1)\times n\) 행렬로 나타낼 수 있다. 그리고 블록 대각 행렬 형태로 식 (15)처럼 각 위성군의 D 행렬을 쌓아 모든 위성군에 대해 변환행렬을 구성한다.

\(D=blkdiag \left [ D^{GPS}, ~ D^{GAL}, ~ D^{BDS} \right ]\)                                                                   (15)

\(D^{GPS}=\left [ \begin{matrix} 1&-1&0 &⋯& 0 \\ 1&0&-1&⋯&0 \\ &⋮& &⋱&⋮\\ 1&0&0&⋯&-1 \end{matrix} \right ]_{(n-1)\times n}\)                                                      (16)

여기서 \(D^{GPS},D^{GAL},D^{BDS}\)는 각 위성군의 이중 차분 변환 행렬이며 위성군의 위성 개수와 pivot 위성에 따라 달라진다.

코드와 반송파 위상의 단일 차분 측정 오차 공분산 행렬 \(R_{\Phi_{sd,w}}\)과 \(R_{\rho_{sd,n}}\)은 D 행렬의 형태와 마찬가지로 다중 위성군 행렬 형태로 구성한다. 이때 코드와 위상 간, L1과 L2 주파수 간의 교차상관관계 (cross-correlation)는 없다고 가정한다.

\(R_{\Phi_{sd,w}}=blkdiag\left[R_{\Phi_{sd,w}}^{GPS}, R_{\Phi_{sd,w}}^{GAL}, R_{\Phi_{sd,w}}^{BDS}\right]\)                                                                (17)

\(R_{\rho_{sd,n}}= blkdiag\left[R_{\rho_{sd,n}}^{GPS}, R_{\rho_{sd,n}}^{GAL}, R_{\rho_{sd,n}}^{BDS}\right]\)                                                                  (18)

\(R_{\rho_{sd,n}}^{GPS}= \left[ \begin{matrix} \left(\sigma_{r,\rho_n}^{1,GPS}\right)^2+\left(\sigma_{b,\rho_n}^{1,GPS} \right)^2&&\\&⋱&\\&&\left(\sigma_{r,\rho_n}^{n,GPS}\right)^2+\left(\sigma_{b,\rho_n}^{n,GPS}\right)^2 \end{matrix} \right]_{n×n}\)                         (19)

\(R_{\Phi_{sd,w}}^{GPS}=\left[ \begin{matrix} \left(\sigma_{r,\Phi_w}^{1,GPS}\right)^2+\left(\sigma_{b,\Phi_w}^{1,GPS}\right)^2&&\\&⋱&\\&&\left(\sigma_{r,\Phi_w}^{n,GPS}\right)^2+\left(\sigma_{b,\Phi_w}^{n,GPS}\right)^2 \end{matrix} \right]_{n×n}\)                         (20)

여기서 \(\sigma_{r,\Phi_w}^s\)와 \(\sigma_{b,\Phi_w}^s\)은 각각 이동국과 기준국의 차분 전 위상 확장파장 조합 측정 오차의 표준편차이며, \(\sigma_{r,\rho_n}^s\)와 \(\sigma_{b,\rho_n}^s\)은 차분 전 코드 축소파장 조합 측정 오차의 표준편차이다. 이는 식 (21), (22)와 같이 계산된다.

\(\sigma_{\Phi_n }=\sqrt{\left({\lambda_{L2} \over \lambda_{L2}-\lambda_{L1} }\right)^2 \sigma_{\Phi_{L1}}^2+\left({\lambda_{L1} \over \lambda_{L2}-\lambda_{L1} }\right)^2 \sigma_{\Phi_{L2}}^2 }\)                                                                   (21)

\(\sigma_{\rho_n }=\sqrt{\left({\lambda_{L2} \over \lambda_{L2}+\lambda_{L1} }\right)^2 \sigma_{\rho_{L1}}^2+\left({\lambda_{L1} \over \lambda_{L2}+\lambda_{L1} }\right)^2 \sigma_{\rho_{L2}}^2 }\)                                                                     (22)

L1과 L2에서의 반송파 위상 측정 오차의 표준편차 \(\sigma_\rho\)와 코드 측정 오차의 표준편차 \(\sigma_\Phi\)를 동일하다 가정하고 정량화된 값을 사용한다.

\(\sigma_{\rho_{L1} }=\sigma_{\rho_{L2} }=\sigma_\rho={0.3 \over \sin⁡(elev)}\)                                                                                     (23)

\(\sigma_{\Phi_{L1} }=\sigma_{\Phi_{L2} }=\sigma_\Phi={0.003 \over \sin⁡(elev)}\)                                                                                   (24)

여기서 elev는 이동국 혹은 기지로부터 측정한 각 위성의 고도각으로, 가중치로 사용된다.

2.2 확장파장 조합 미지정수 결정

이동국과 기준국의 정밀한 상대 벡터 측위는 반송파 위상 이중 차분 측정치의 신뢰 높은 미지정수의 추정이 필수적이다. 앞서 설명한 확장파장 조합과 축소파장 조합 (식 (9), (10))으로부터 이중 차분을 하고 이를 선형화하여 상대 벡터 형태 (식 (25), (26))를 구성하고, 반송파 위상 이중 차분 미지정수의 해를 LAMBDA 기법을 통해 계산한다.

\(\Phi_{dd_w}=Gb+ \Lambda_w N_{dd,w}+\varepsilon_{\Phi_{dd,w}}\)                                                                       (25)

\(\rho_{dd_n}=Gb+\varepsilon_{\rho_{dd,n}}\)                                                                                                (26)

여기서 \(G\)는 각 위성들과 pivot 위성 간 정규화된 시선각 벡터의 차분 행렬이며, \(b\)는 기준국에 대한 이동국의 상대 위치 벡터, \(\Lambda_w\)는 확장파장 조합의 파장 블록 대각 행렬이고, \(N_{dd,w}\)는 이중 차분 확장파장 조합의 미지정수 벡터이다.

LAMBDA 기법으로 미지정수를 결정하기 위해 최소제곱법 (Least-squares)의 실수해 (floating solution)와 공분산 행렬을 이용한다. \(\hat{b}\)와 \(\hat{N}\)은 단일 epoch에서 식 (27-30)과 같이 계산된다.

\(\hat{Q}_b=\left(G^T R_{\rho_{dd,n}}^{-1} G\right)^{-1}\)                                                                                     (27)

\(\hat{b}=\hat{Q}_b G^T R_{\rho_{dd,n}}^{-1} \rho_{dd,n}\)                                                                                  (28)

\(\hat{Q}_N=\Lambda_w^{-1} \left(R_{\Phi_{dd,w}}+G\hat{Q}_b G^T\right)\Lambda_w^{-1}\)                                                        (29)

\(\hat{N}_{dd,w}=\Lambda_w^{-1} \left(\Phi_{dd,w}-G\hat{b}\right)\)                                                                       (30)

여기서 \(\hat{Q}_b\)와 \(\hat{Q}_N\)는 각 추정치의 공분산 행렬이다. LAMBDA에서 나온 미지정수 \(\check{N}\)가 일정 수준 이상의 신뢰도를 갖는 경우 미지정수가 결정 (fixed)되며, 상대벡터 \(\check{b}\)와 공분산 행렬인 \(\check{Q}_b\)는 식 (31), (32)와 같이 계산된다.

\(\check{Q}_b=\left(G^T \left(R_{\rho_{dd,n}}^{-1}+R_{\Phi_{dd,w}}^{-1}\right )G\right)^{-1}\)                                                                          (31)

\(\check{b}=\check{Q}_b G^T \left(R_{\rho_{dd,n}}^{-1} \rho_{dd,n}+R_{\Phi_{dd,w}}^{-1} \left(\Phi_{dd,w}-\Lambda\check{N}\right)\right)\)                                           (32)

최소제곱법 대신 칼만필터가 사용된 경우 수신기의 위치, 속도, 단일 차분 미지정수를 상태벡터로 가지며 식 (33-37)을 통해 계산된다.

\(\hat{x}_{k+1} (t_i )=F_k^{k+1} \hat{x}_k (t_{i+1} )\)                                                                                                              (33)

\(P_{k+1} (t_i)=F_k^{k+1} P_k (t_{i+1}) F_k^{k+1T}+Q_k^{k+1}\)                                                                               (34)

\(\hat{x}_k (t_{i+1} )=\hat{x}_k (t_i )+K_k (y_k-h(\hat{x}_k (t_i )))\)                                                                               (35)

\(P_k (t_{i+1} )=(I-K_k H(\hat{x}_k (t_i )))P_k (t_i )\)                                                                                      (36)

\(K_k=P_k (t_i )H(\hat{x}_k (t_i )) \left(H(\hat{x}_k (t_i )) P_k (t_i )H(\hat{x}_k (t_i ))^T+R_k \right )^{-1}\)                           (37)

여기서 \(\hat{x}_k, P_k\)는 각각 상태벡터와 상태벡터에 해당하는 공분산행렬이고, \(h(x), H(x), R_k\)는 측정치 모델, 측정치 오차 모델을 나타낸다. \(F_k^{k+1}\)와 \(Q_k^{k+1}\)는 각각 상태 전이행렬과 시스템 잡음 공분산행렬이다. 칼만필터를 통해 계산된 상태벡터와 공분산행렬로부터 미지정수 float \(( \hat{x}'_k )\)과 공분산 행렬\((Q_N)\)은 식 (38-40)와 같이 계산된다.

\(\hat{x}'_k=G\hat{x}_k (t_{i+1} )=\left(\hat{r}_r^T, \hat{v}_r^T, \hat{N}^T \right)^T\)                                                                             (38)

\(P'_k=GP_k (t_{i+1} ) G^T=\left[ \begin{matrix} Q_R&Q_{NR}\\Q_{RN}&Q_N \end{matrix} \right]\)                                                                     (39)

\(G=\left[ \begin{matrix} I_{6×6}&\\&D\end{matrix} \right]\)                                                                                                               (40)

2.3 Moving Baseline RTK 구현을 위한 측정값 구성 및 상태변수 지정

2.1절과 2.2절에 나열된 수식들은 고정된 기준국 수신기를 이용한 RTK 기법과 유사하나 원활한 moving baseline RTK 구현할 때 주의해야 할 점은 다음과 같다.

\(\cdot\) 측정값의 시간 동기화: 움직이는 수신기들의 측정값으로 확장파장 조합과 이중 차분 조합을 구성할 때 측정값들이 시간동기화가 되어있는지 확인하여야 한다. RTK 기능이 내재되어 있는 고가의 수신기들은 이미 지정된 시간을 기준으로 각 수신기들이 원시 측정값들을 생성하여 시간 동기화가 필요 없지만 저가의 수신기를 사용할 때는 측정값의 time-of-week의 차이를 확인하고 측정값들의 시간동기화를 진행하여야 한다.

\(\cdot\) 상태변수지정: 고정된 기준국의 정밀한 위치를 사용하는 RTK 기법에서 이동국의 절대 위치를 변수로 지정하는 경우가 있다. Moving baseline RTK 문제에서 이동국의 절대 위치를 변수로 사용할 경우 (특히 칼만 필터에서) 기준국의 코드 기반 절대위치를 사용해야 하며 이 때 원치 않는 오차가 추가로 필터에 유입될 수 있다. 이를 방지하기 위해 상대벡터를 변수로 지정하는 것이 유리하다.

\(\cdot\) Fixed baseline 해 측정: LAMBDA에서 도출된 미지정수해가 ratio test를 통과할 경우 고정해 (fixed solution)라 판단을 하며, 이를 이용하여 고정해 기선 (fixed baseline)을 계산한다. 고정해 기선계산 시 미지정수해가 정확히 도출되었다 하더라도 실수해 기선 (float baseline) 정확도가 고정해 기선에 영향을 미칠 수 있다. Moving baseline RTK 경우 움직이는 두 수신기의 가속도가 급격히 바뀔 때 칼만필터가 이를 적절히 추정하지 못할 수 있으며 이러한 오차는 실수해 기선에 포함되어 결국 고정해 기선의 오차도 따라 증가할 수 있다. 따라서 이전 위치측정 결과에 독립적인 고정해 기선을 결정하는 식 (32)를 사용하는 것이 바람직하다.

3. 실험 환경 및 구성

3.1 ZEP-F9P

본 실험에서 사용한 수신기는 u-blox의 ZED-F9P 모듈이며, ZED-F9P는 다중 위성군과 다중 주파수를 지원한다. 지원 가능한 위성군은 GPS, GLONASS, Galileo, BeiDou, Quasi-Zenith Satellite System (QZSS)이 있으며, 총 184개 채널을 지원한다. 측위 정확도는 0.01 m + 1 ppm Circular Error Probability (CEP)이며, 최대 20 Hz까지 위치측정이 가능하고, 수신기 크기는 110 x 55 mm이다. 안테나는 이중 주파수 신호 수신이 가능한 ANN-MB series GNSS Antenna를 사용하였다. Fig. 1은 본 연구에서 알고리즘 평가를 위해 사용한 안테나 및 수신기다.

Fig. 1. GNSS antenna and receiver: ANN-MB series & ZED-F9P.

3.2 실험 환경

본 연구는 기준국과 이동국이 함께 군집 비행하는 moving baseline 실험을 위해 2대의 드론을 동시에 띄워 실험을 진행하였다. 실험은 홍익대학교 화성캠퍼스 운동장에서 총 두 번 실시하였으며, DJI 사의 F450에 ZED-F9P를 탑재하여 GNSS 데이터를 수신하였다. 첫번째 실험은 약 5분, 두번째 실험은 약 12분간 진행하였다. 데이터 수신 시 이동국과 기준국의 데이터를 각각 수신하였으며, 모든 데이터는 1초 간격으로 수신하였다. Fig. 2는 군집 비행을 실험한 홍익대학교 화성 캠퍼스 운동장이다. Fig. 3은 본 실험에서 사용한 GNSS 데이터 수집 및 군집 비행용 드론이다.

Fig. 2. Snapshot of Hongik University field & flight path of Drones.

Fig. 3. F9P GNSS receiver and Raspberry Pi mounted on Drone.

4. 상대항법 결과 비교

Fig. 4는 각각 Test 1과 2에 대한 sky plot이다. L1, L2 신호 모두 수신되는 동시에 cut-off 각도인 25도가 넘는 위성을 표시하였고, 해당 위성들을 RTK 측위에 사용하였다. 25도 이하의 위성을 포함하여 RTK 측위에 사용했을 때 residual의 오차가 커지거나 fix가 되지 않는 구간의 길이가 증가하여 cut-off 각도를 25도로 설정하였다. Table 1은 Test 1과 2에서 전체 비행시간 동안 이동국의 평균 GDOP를 각 위성군에 따라 결과를 나타낸 것이다. 위성이 제일 많은 test2의 GPS/GAL/BDS 위성군에서 평균 GDOP가 1.24로 제일 낮고 위성수가 감소할수록 GDOP가 높아지므로 모든 위성군을 사용한 경우 각 실험 내에서 가장 좋은 위성환경을 가진다.

Fig. 4. Skyplot of satellites: (left) Test 1, (right) Test 2.

Table 1. Average GDOP of each satellite constellation combinations.

Satellite constellations	Test 1		Test 2
	Average GDOP	# of SV	Average GDOP	# of SV
GPS/GAL/BDS GPS/GAL GPS/BDS	1.64 1.79 3.14	9 6 5	1.24 1.45 2.16	10 7 6

 

또한 확장파장 조합과 이중 주파수 측정값을 사용하여 5가지 RTK 측위 성능을 Tables 2-7에 나타냈다. Tables 2-4은 Test 1의 결과이며, Table 5-7은 Test 2의 결과이다. Tables 2과 5는 GPS/GAL/BDS 위성군을 모두 사용한 경우이며, Tables 3과 6는 GPS/GAL, Tables 4와 7은 GPS/BDS 위성군 조합 결과를 나타낸다. Figs. 5-9는 Test 2의 모든 위성군을 사용한 RTK 측위 기법에 따른 LAMBDA의 ratio test 결과 및 고정해와 실수해로 계산된 이동국과 기지국의 상대벡터 오차를 나타낸다. 이때 ratio test의 임계값을 3으로 설정하였으며, ratio test를 통과한 미지정수해를 식 (32)에 적용하여 고정해를 계산하였다. Tables 2-7에 제시된 성능지표는 TTFF, 미지정수해 도출 성공 확률, 위치오차의 RMS 값이다. 표에서 fixed ratio는 전체 실험 시간 중 LAMBDA ratio를 통과한 비율이며, success ratio는 전체 고정해 중 참 값과 비교했을 때 올바른 값을 가진 비율 값이다. 상대 벡터와 미지정수해 도출의 참 값으로는 오픈소스 RTK 소프트웨어인 RTKLIB에서 계산된 실험 지역 주변 고정된 기준국과의 RTK 결과를 사용하였다 (Takasu & Yasuda 2009). 반올림 기법을 사용한 확장파장 조합은 최소 50 epoch의 실수 미지정수해를 batch로 반올림하여 미지정수해를 도출하였으며, 50 epoch 이전은 실수 미지정수해를 사용하여 상대 벡터를 계산하였다. 50 epoch을 기준으로 반올림한 이유는 해당 epoch일 때 본 실험에서 모든 위성에 대한 확장파장 조합의 잡음 불확도 (noise uncertainty) 가 0.1 cycle 미만으로 낮아졌기 때문이다.

Table 2. Moving baseline RTK results with GPS / GAL /BDS (Test 1).

GPS / GAL / BDS			Widelane KF	Widelane LS	Widelane round-off	Dual freq. KF	Dual freq. LS
# of SV used TTFF Fixed ratio Success ratio			9 1 1 1	9 1 0.60 1	9 50	9 1 1 1	9 1 1 1
RMSE (cm)	Fix	E N U All	5.15 3.15 21.46 22.30	4.70 2.89 20.08 20.83	5.52 3.30 22.07 22.99	0.78 0.59 4.72 4.83	0.77 0.52 4.71 4.80
	Float	E N U All	NA	68.99 92.44 409.77 425.70	20.73 31.56 90.45 98.01	NA

 

Table 3. Moving baseline RTK results with GPS / GAL (Test 1).

GPS / GAL			Widelane KF	Widelane LS	Widelane round-off	Dual freq. KF	Dual freq. LS
# of SV used TTFF Fixed ratio Success ratio			6 2 0.89 1	6 2 0.20 0.74	6 50	6 1 1 1	6 1 0.52 1
RMSE (cm)	Fix	E N U All	6.45 3.42 24.68 25.73	24.52 53.38 77.90 97.56	6.57 3.51 24.79 25.88	0.81 0.54 4.46 4.57	0.75 0.44 4.69 4.77
	Float	E N U All	39.40 55.22 88.71 111.67	71.24 88.81 359.87 377.45	38.34 71.34 171.16 189.34	NA	73.26 98.79 371.48 391.31

 

Table 4. Moving baseline RTK results with GPS /BDS (Test 1).

GPS / BDS			Widelane KF	Widelane LS	Widelane round-off	Dual freq. KF	Dual freq. LS
# of SV used TTFF Fixed ratio Success ratio			5	5	5	5	5
RMSE (cm)	Fix	E N U All	Lack of SV
	Float	E N U All

 

Table 5. Moving baseline RTK results with GPS / GAL /BDS (Test 2).

GPS / GAL / BDS			Widelane KF	Widelane LS	Widelane round-off	Dual freq. KF	Dual freq. LS
# of SV used TTFF Fixed ratio Success ratio			10 1 1 1	10 1 0.73 1	10 50	10 1 1 1	10 1 0.95 1
RMSE (cm)	Fix	E N U All	3.96 4.15 13.55 14.72	3.90 4.08 13.18 14.34	5.58 12.86 26.58 30.05	0.49 0.48 1.46 1.62	0.48 0.43 1.43 1.57
	Float	E N U All	NA	58.39 53.29 262.31 273.96	9.15 36.73 121.61 127.36	NA	62.39 47.14 317.74 327.22

 

Table 6. Moving baseline RTK results with GPS / GAL (Test 2).

GPS / GAL			Widelane KF	Widelane LS	Widelane round-off	Dual freq. KF	Dual freq. LS
# of SV used TTFF Fixed ratio Success ratio			7 1 0.93 1	7 1 0.39 0.95	7 50	7 1 0.99 1	7 1 0.54 1
RMSE (cm)	Fix	E N U All	5.59 4.84 16.13 17.74	19.04 19.99 106.0 109.53	8.72 14.75 19.94 26.29	0.52 0.44 1.68 1.81	0.51 0.40 1.60 1.73
	Float	E N U All	8.58 47.22 120.05 129.28	61.53 53.73 245.65 258.87	19.20 53.54 130.99 142.80	12.27 38.60 136.59 142.47	65.24 51.19 255.34 268.46

 

Table 7. Moving baseline RTK results with GPS / BDS (Test 2).

GPS / BDS			Widelane KF	Widelane LS	Widelane round-off	Dual freq. KF	Dual freq. LS
# of SV used TTFF Fixed ratio Success ratio			6 2 0.97 1	6 2 0.40 0.96	6 50	6 2 0.99 1	6 1 0.46 1
RMSE (cm)	Fix	E N U All	13.53 7.62 55.89 58.00	48.59 25.06 168.25 176.91	13.67 7.76 56.44 58.59	2.20 1.40 8.28 8.69	2.29 1.49 7.97 8.42
	Float	E N U All	67.55 10.59 505.91 510.51	106.82 63.01 441.01 458.12	60.84 24.61 412.84 418.02	20.98 5.95 112.34 114.44	94.96 54.48 379.35 394.83

 

Fig. 5. Relative position errors of (left) rounded solution and (right) float solution of widelane round-off method using GPS/GAL/BDS satellite constellation for Test 2.

Fig. 6. Relative position errors of (left) fixed solution and (middle) float solution of widelane Kalman filter method using GPS/GAL/BDS satellite constellation for Test 2, and (right) ratio test of this method.

Fig. 7. Relative position errors of (left) fixed solution and (middle) float solution of widelane least square method using GPS/GAL/BDS satellite constellation for Test 2, and (right) ratio test of this method.

Fig. 8. Relative position errors of (left) fixed solution and (middle) float solution of dual frequency least square method using GPS/GAL/BDS satellite constellation for Test 2, and (right) ratio test of this method.

Fig. 9. Relative position errors of (left) fixed solution and (middle) float solution of dual frequency Kalman filter method using GPS/GAL/BDS satellite constellation for Test 2, and (right) ratio test of this method.

확장파장 조합이 널리 사용되는 이유는 확장된 파장 길이로 미지정수해가 보다 안정적으로 도출될 수 있다고 여겨지기 때문이지만 Teunissen (1997)에서 증명된 바와 같이 확장파장 조합은 LAMBDA의 측정값 decorrelation 프로세스의 한 과정에 불과하며 이중 주파수 측정치를 사용한 LAMBDA 기법이 확장파장 조합 보다 미지정수해 도출 성공률이 더 높다. 이러한 특성은 Tables 3, 6, 7의 확장파장 조합이 이중 주파수 보다 미지정수해 도출 성공 확률이 낮게 나타나는 것에서도 확인된다. 이중 주파수의 도출 성공 확률이 모두 1, 즉 정확히 도출되었지만 확장파장 조합의 최소제곱법의 도출 성공 확률이 최소 74%이다 (Table 3). Tables 2, 5와 Fig. 6에서와 같이 확장파장 조합을 사용하였을 때 위치측정 오차가 동일 위성군의 이중주파수 RTK 보다 5배 이상 높으며, 이는 확장파장 조합에서 파장의 길이가 증가하여 노이즈가 증폭되었기 때문이다. 또한 확장파장 조합은 이중 주파수 측정값이 반드시 요구되므로 이러한 조건이 수반되지 않는 최소제곱법과 칼만필터를 이용한 moving baseline RTK보다 가용위성 개수가 적을 수 있다. 따라서 이중 주파수 기반 최소제곱법이나 칼만필터를 사용하는 것이 미지정수해 도출과 상대위치벡터 산출에 있어 바람직하다 할 수 있으나, 이 두가지 방법은 LAMBDA 기법에서 ratio test의 임계값 (threshold)을 넘거나 부트스트래핑 (boot-strapping) 기반의 미지정수해 도출 성공확률 임계값을 통과해야만 미지정수해 도출 결과를 신뢰할 수 있고, 정밀한 상대위치 벡터를 산출할 수 있다. 본 논문에서는 ratio test를 사용하였으며 임계값은 3으로 설정하였다. 가용 위성군 또는 가용위성 개수가 적을 경우 임계값을 통과하지 못할 수 있으며 이러한 경우 실수 미지정수해 기반 상대위치벡터를 사용해야 하기 때문에 위치측정 정확도의 저하를 감수해야 한다. 확장파장 조합이 LAMBDA보다 미지정수해 도출 성공확률은 우수하지 못하나 LAMBDA 솔루션이 threshold를 통과하지 못할 경우에 임시적으로 상대벡터 도출에 사용될 수 있고 실수해 기반 상대위치 벡터 보다 정밀하다 (Misra & Enge 2006).

Tables 2-7에서 확장파장 조합 칼만필터 기법의 fixed ratio는 약 0.9 이상, 이중주파수의 칼만필터 또한 1에 근접한 수치를 보이는 반면 최소제곱법의 fixed ratio는 위성군을 모두 사용했을 때 (Tables 2, 5) 0.7 이상이며, 나머지 두 조합의 위성군 결과 (Tables 3, 4, 6, 7)에서는 앞선 결과보다 낮은 0.2~0.5이다. 따라서 가용 위성의 개수가 높을수록 LAMBDA에서 미지정수가 풀릴 가능성이 높다. Test 1의 GPS/BDS 위성군 조합의 RTK 결과인 Table 4에서는 가용 위성 수가 5개로, LAMBDA에서 미지정수가 풀리지 않는 경우가 급증하였으며 미지정수가 풀리더라도 RTK 해의 오차가 수 km에 이를 정도로 매우 부정확한 결과를 보여 표에는 반영하지 않았다. 해당 위성군의 평균 GDOP는 3.14로 동일 실험의 GDOP 중 가장 높은 수치로 측정치 정밀도가 가장 낮다고 볼 수 있다 (Table 1). 반면 동일 실험의 GPS/GAL 위성군을 사용한 결과 (Table 3)에서는 GPS/BDS 위성군보다 가용 위성이 하나 더 많은 6개지만 모든 RTK 기법에 대해 항법해를 도출 가능하다.

또한 칼만필터를 사용했을 때 모든 moving baseline RTK결과 (Tables 2~7)에서 최소제곱법보다 높은 fixed ratio와 성공률을 보였으며 칼만필터 방식의 RTK 기법이 위성수가 6개인 환경에서 (Table 7) 최소제곱법 보다 두 배 가량 고정해를 도출하였다. 이를 통해 GDOP가 상대적으로 높은 위성환경에서도 칼만필터가 미지정수해 도출의 이점을 가진다고 볼 수 있으나, GDOP의 상승으로 인해 측정치 노이즈가 증가하여 고정 미지정수해의 상대벡터 오차가 커진다.

GPS/GAL/BDS 위성군을 모두 사용한 Tables 2, 5에서는 이중 주파수 최소제곱법 조합의 fix 미지정수해 기반 상대위치 벡터가 2~4 cm 가량 위치측정 오차를 가지지만 float 미지정수해의 위치 오차의 경우 fix일 때 보다 40배 이상 가까이 오차가 증가하게 된다. 반면 동일 위성군의 widelane의 round-off 조합의 상대위치 벡터의 위치 측정 결과는 평균적으로 약 30 cm 정도의 오차가 존재하지만 동일 위성 비교군의 다른 float 미지정수해 기반 상대위치 벡터 오차의 RMS보다는 작음을 확인할 수 있다. Fig. 10은 Fig. 8의 결과에서 이중 주파수 최소제곱법의 미지정수해가 float이 되는 해당 epoch에 실수 미지정수해 대신 확장파장 조합의 반올림 기법을 통해 도출한 미지정수해로 계산한 상대 벡터의 위치 측정 오차이다. Table 8은 해당 epoch에서의 이중주파수 최소제곱법 실수해와 확장파장 조합 반올림 기법의 미지정수해 RMSE의 결과이다. 반올림 기법의 미지정수를 사용했을 때 실수 미지정수해의 위치 측정 오차보다 평균 1/10 가량 감소함을 확인할 수 있으며, 실수 미지정수해의 대체 항법해로서 사용이 유의미함을 보인다.

Fig. 10. Comparison of relative position errors of float solution of dual frequency least square using GPS/GAL/BDS satellite constellation and widelane rounded solution at the same epoch in Test 2.

Table 8. RMSE of widelane round-off and dual frequency least square with GPS/GAL/BDS at same epochs (Test 2).

RMSE (cm)	Widelane round-off	Dual freq. least squares
E N U All	5.55 8.58 27.56 29.39	62.39 47.14 317.74 327.22

 

5. 결론

본 연구에서는 이동국과 기준국이 함께 이동하는 moving baseline 상황에서 확장파장 조합과 이중 주파수를 사용하였을 경우에 대해 반올림, LAMBDA, 칼만 필터 기법을 이용한 RTK 기법들의 결과를 비교하였다. 앞선 상대항법 결과 비교에서는 평균 GDOP가 가장 낮은 GPS/GAL/BDS 위성군을 다른 비교군 보다 위성환경이 좋음을 보였으며, GDOP가 낮고 위성개수가 높아질수록 Fixed ratio와 success ratio는 높으며, rmse는 감소하여 측위 성능이 향상된다. 또한 측정치 조합 측면에서는 Test 1과 2 모두 확장파장 조합 RTK 기법의 fixed ratio가 이중 주파수에서 동일한 기법의 fixed ratio보다 대부분 0.06~0.40 가량 낮고, 고정해의 RMS 오차 또한 10 cm 이상 크기 때문에 이중주파수 측정치 조합이 확장파장 측정치 조합보다 RTK 성능이 좋다고 볼 수 있다. Tables 2-7의 결과 모두 칼만필터가 동일 측정치 조합의 최소제곱법에 비해 상대벡터 오차, 미지정수 고정 비율, 도출 성공확률 모두 향상된 성능을 보이며, 위성 수가 적어 상대적으로 GDOP가 높은 Table 7에서 fix RMSE가 1 m 이상 차이나는 등 GDOP가 클수록 차이가 크다. 따라서 가시 위성 수가 최소 6 이상이며 GDOP가 2 이하의 위성환경에서 moving baseline의 정밀한 상대벡터 값의 계산을 위해 이중 주파수 최소제곱법 조합이나 칼만 필터 조합을 쓰는 것이 합리적이다. 하지만 LAMBDA에서 도출된 미지정수해의 결과가 ratio test를 넘지 못하는 경우에 한해서는 확장파장 조합의 반올림 기법의 상대 위치 벡터 측정 오차가 이중 주파수의 실수해 기반의 측정 오차보다 작기 때문에 이중 주파수 RTK 기법과 확장파장 조합의 반올림 기법이 동시에 사용되는 것을 권장한다. 이 경우 두 개의 RTK 기법을 동시에 운용함에 따른 연산량의 부담을 감당할 수 있는 프로세서가 요구된다.