GPS/GF-INS Integrated Navigation System with High Rate Position, Velocity, and Attitude Aiding of GPS

Abstract

In this paper, a GPS/GF-INS integrated navigation system is proposed, in which the high rate attitude aiding signal, the high rate position and velocity aiding of GPS receiver is used for the cube structure of the GF-IMU, effectiveness of the proposed GPS/GF-INS integrated navigation system was shown when the vehicle follows two trajectories, circling and spiraling. Performance evaluation results show that the proposed GPS/GF-INS integrated navigation method gives better navigation outputs when the attitude output of GPS is used and more better navigation outputs are obtained when the rate of GPS aiding signal is higher.

1. INTRODUCTION

항법은 항체의 위치, 속도, 자세를 구하는 것이며, 관성항법 시스템 (Inertial Navigation System, INS)은 가속도계 3개와 자이로 3개로 이루어진 Inertial Measurement Unit (IMU)의 측정치로부터 항체의 위치, 속도, 자세를 제공하는 시스템이다 (Titterton & Weston 2004). 그런데, 6개의 가속도계만을 이용해도, 종래의 IMU와 같이 3차원의 가속도와 3차원의 각가속도를 측정할 수 있다 (Chen et al. 1994). 따라서, 항체가 자이로스코프의 측정범위를 넘는 고속 회전을 할 경우나 복잡하고, 가격이 비싸며, 전력을 많이 소비하는 자이로스코프의 사용을 피하기 위하여 가속도계만으로 항법 시스템을 구성할 수 있다 (Tan & Park 2005, Pachter et al. 2013). 항법급 IMU인 HG9900에 사용되는 GG1320 자이로는 1.6 W를 소모하며, QA2000 가속도계는 480 mW를 소모한다 (Honeywell 2004, 2021). 일반적으로 직선운동과 회전운동을 하는 항체의 1축의 각속도를 구하려면, 가속도계 2개가 필요하므로, 자이로가 약 1.7배의 전력을 소모하는 것을 알 수 있다. 또한, 원점으로부터 일정 거리만큼 떨어진 지점에 가속도계를 배치하여, 각가속도와 각속도를 얻을 수 있는데, 각속도 추정 성능은 가속도계의 원점으로부터의 거리에 따라 달라진다 (Pachter et al. 2013). 즉, 거리를 키울 경우, 각속도 추정 시 가속도계 오차에 대한 영향이 줄어들기 때문에, 동급의 자이로보다 정확한 각속도를 얻을 수 있다. 가속도계만으로 항법시스템을 구성할 때, 항체에 배치한 가속도계를 Gyro Free-IMU (GF-IMU)라 하고, GF-IMU의 출력을 이용하여 구성한 항법 시스템을 Gyro Free- INS (GF-INS)라고 한다.

Fig. 1. GF-INS.

GF-IMU를 사용하는 개념은 1975년도에 제안되었으며, 그동안은 MEMS 기술을 활용하는 방안으로 제안되었다 (Qin et al. 2009, Pachter et al. 2013). 즉, MEMS 가속도계가 MEMS 자이로스코프에 비하여 더 소형이며, 가격이 싸고, 전력을 덜 소비하는 장점 때문에, GF-IMU를 유도무기, 무인 항공기에 이용하려는 연구를 시도하였다 (Liu et al. 2019). 그런데, MEMS 가속도계는 종래의 광 기술을 이용한 가속도계나 항법급 가속도계에 비하여 성능이 매우 나쁘다. 더욱이 GF-INS는 GF-IMU의 출력으로부터 각가속도를 계산하고, 이로부터 각속도를 구하므로, 종래의 IMU로부터 구한 INS의 항법 결과에 비하여 오차의 증가 속도가 훨씬 빠르다. 최근에는 냉각 원자 간섭에 의하여 항체의 가속도를 매우 정확히 측정하는 것이 가능하면, Global Positioning System (GPS)의 도움 없이 신뢰성 있는 항법 결과를 얻을 수 있음을 예상하고 있다 (McGuirk et al. 2002, Sorrentino et al. 2010, Pachter et al. 2013, Nusbaum & Klein 2017). GF-INS는 Fig. 1에서 나타낸 바와 같이 각가속도와 비력을 구하여 각속도를 추정하고, 각속도와 가속도를 적분하여 항법 결과를 얻으므로, INS와 마찬가지로 센서, 즉, 가속도계의 오차에 의하여 오차가 시간에 따라 증가한다. 그리고, INS에 비하여 각가속도를 적분하는 과정이 더 있으므로, 항법 오차는 더욱 더 빨리 증가한다.

GF-INS의 항법 오차를 줄이기 위해 보조 신호 (aiding signal)로 GPS의 출력을 이용하는 GPS/GF-INS 통합항법 시스템을 생각할 수 있다. Park & Tan (2002)은 정육면체 각 면에 1개씩의 가속도계를 배치한 GF-IMU에 대하여 GF-INS의 선형과 비선형 오차모델을 도출하고, 약결합 방식의 GPS/GF-INS 통합항법 시스템을 구성하였다. 그리고, GF-INS 단독항법과 GPS/GF-INS 통합항법 시스템의 성능을 비교하였다. 선형과 비선형 모델 각각을 이용한 2가지의 통합 칼만필터를 구성하였으며, 특정한 궤적에 대하여 GPS/GF-INS 통합항법 시스템 위치 결과가 GF-INS 단독 위치 결과에 비하여 좋은 결과를 제공함을 보였다. Marques Filho et al. (2006)은 Park & Tan (2002)과 같은 GF-IMU 구성에 대하여 유사한 방식의 약결합 GPS/GF-INS 통합항법 시스템의 시뮬레이션을 수행한 결과, 통합항법 시스템의 결과가 GF-INS 항법 결과에 비하여 우수하며, 항법 오차가 GPS 측정치의 표준편차 (1 σ) 이내임을 보였다. Son et al. (2020)도 동일한 구조로 GPS/GF-INS 통합항법 시스템을 구성하였으며, 가속도계 오차가 잡음뿐인 경우, 통합항법 시스템의 위치, 속도, 자세 오차가 발산하지 않음을 보였다. Edwan et al. (2009)은 원점과 각 축에 가속도계 triad를 배치한 12개의 가속도계로 이루어진 GF-IMU에 대하여, GPS/GF-INS 통합항법 시스템을 구성하고, 통합 칼만필터의 측정치가 위치, 속도 차이인 경우, 위치, 속도, 각속도 차이인 경우에 대해 위치와 자세 추정 성능을 비교하였다. 각속도 측정치를 추가한 경우에는 위치 추정 성능이 크게 향상되며, 자세 오차의 발산 속도를 늦출 수 있음을 보였다. 그런데, 전술한 결과들을 보면, 가속도계 오차의 가정이 실제와 거리가 있으며, 위치, 속도, 자세의 모든 항법 결과가 수렴하는 결과들은 제시하지 않고 있다. 즉, 일반적인 GPS의 위치, 속도를 보조 신호로 이용해서는 GPS/GF-INS 통합항법의 원하는 결과를 얻을 수 없음을 알 수 있다.

본 논문에서는 위치, 속도뿐만 아니라 자세 정보까지 제공하는 GPS 수신기와 보조 신호의 출력을 높이는 방법을 이용하는 새로운 GPS/GF-INS 통합항법 시스템을 제안하고, 그 성능을 제시하고자 한다. 앞으로 2절에서는 GPS/GF-INS 통합항법 알고리즘을 설계하기 위해 필요한 GF-INS의 항법방정식, GF-INS 오차모델, 통합 칼만필터에 대해 상세히 서술한다. 3절에서는 GPS 수신기의 출력율과 위치, 속도 측정치에 자세 측정치를 추가하여 확인한 성능평가 결과를 제시한다. 마지막으로 4절에서 본 연구의 결론을 정리하고, 추후 계획을 서술한다.

2. 위치, 속도 및 자세 정보를 제공하는 GPS와 GF-INS의 통합

본 논문에서 새롭게 제안하는 GPS/GF-INS 통합항법 시스템은 Fig. 2와 같이 위치, 속도, 자세를 제공하는 GPS 수신기, GF-IMU, GF-INS, 통합 칼만필터로 구성된다. GF-INS의 항법 결과는 통합 칼만필터에서 추정한 GF-INS 항법 오차를 이용하여 보정한다.

Fig. 2. Proposed GPS/GF-INS integrated navigation system.

2.1 GF-INS 항법 방정식

Fig. 1에서 보듯이, GF-INS는 가속도계 측정치로부터 계산한 각가속도를 적분하는 각속도 추정기와 INS로 구성된다. 따라서, GF-INS의 항법 방정식은 식 (1-3)에 제시한 INS의 자세, 속도, 위치 미분방정식에 각속도 미분방정식이 추가됨을 알 수 있다 (Klein 2015, Nusbaum & Klein 2017).

\(\mathbf{C} ̇_b^n=\mathbf{C}_b^n \boldsymbol{\Omega}_{ib}^b-(\boldsymbol{\Omega}_{ie}^n+\boldsymbol{\Omega}_{en}^n ) \mathbf{C}_b^n\)                                                                           (1)

\(\dot{\mathbf{v}}^n=\mathbf{C}_b^n \mathbf{f}^b-(2\boldsymbol{\Omega}_{ie}^n+\boldsymbol{\Omega}_{en}^n ) \mathbf{v}^n+\mathbf{g}_l^n\)                                                                 (2)

\(\dot{\mathbf{p}}^n=\left [ \begin{matrix} {\nu_N^n \over \mathbf{R}_M+h} & {\nu_E^n \over (\mathbf{R}_N+h) \cos ⁡\varphi} & -\nu_D^n \end{matrix} \right]^T\)                                                               (3)

여기서 \(\mathbf{C}^n_b\)는 동체좌표계에서 항법좌표계로의 Direction Cosine Matrix (DCM), \(\boldsymbol{\Omega}^n_{ie}\)는 항법좌표계에서 나타낸 지구 자전 각속도, \(\boldsymbol{\Omega}^n_{en}\)는 항법좌표계에서 나타낸 craft rate, \(\mathbf{g}^n_l\)는 항법좌표계에서 나타낸 local gravity, \(\nu_N^n\)는 north축 속도, \(\nu_E^n\)는 east축 속도, \(\nu_D^n\)는 down축 속도, \(\mathbf{R}_M\)은 위도 방향의 반지름 (meridian radius), \(\mathbf{R}_N\)는 경도 방향의 반지름 (normal radius), \(h\)는 고도, \(\varphi\)는 위도를 의미한다.

Fig. 3. N accelerometers in the vehicle.

Fig. 3과 같이 항체에 N개의 가속도계가 부착되어 있다고 할 때 Fig. 3의 i번째 가속도계의 출력은 식 (4)와 같이 나타낼 수 있다 (Schuler et al. 1967, Chen et al. 1994, Park & Tan 2002).

\(y_i=\mathbf{d}_i^{\rm T} \mathbf{f}^b-\mathbf{d}_i^{\rm T} [\mathbf{r}_i^b\times] \boldsymbol{\omega} ̇_{ib}^b+\mathbf{d}_i^{\rm T} (\boldsymbol{\omega}_{ib}^b\times(\boldsymbol{\omega}_{ib}^b\times\mathbf{r}_i^b ))\)                                                   (4)

여기서, \(\mathbf{d}_i\)는 i번째 가속도계의 측정 방향 벡터, \(\mathbf{f}^b\)는 비력, \(\mathbf{r}_i^b\)는 i번째 가속도계의 위치 벡터, \(\boldsymbol{\omega} ̇_{ib}^b\)는 각가속도, \(\boldsymbol{\omega}_{ib}^b\)는 각속도를 의미한다. 식 (4)를 N개에 대해 정리하면, 식 (5)가 된다.

\(\mathbf{Y}=\left[ \begin{matrix} y_1\\ ⋮\\ y_N \end{matrix} \right] =\left[ \begin{matrix} (\mathbf{r}_1^{\rm b}×\mathbf{d}_1 )^{\rm T}&\mathbf{d}_1^{\rm T}\\⋮&⋮\\(\mathbf{r}_N^{\rm b}×\mathbf{d}_N )^{\rm T}&\mathbf{d}_{\rm N}^{\rm T} \end{matrix} \right] \left[ \begin{matrix} \boldsymbol{\omega} ̇_{ib}^b \\ \mathbf{f}^b \end{matrix} \right] +\left[ \begin{matrix} \mathbf{d}_1^{\rm T} (\boldsymbol{\omega}_{ib}^b×(\boldsymbol{\omega}_{ib}^b×\mathbf{r}_1^b )) \\ ⋮ \\ \mathbf{d}_N^{\rm T} (\boldsymbol{\omega}_{ib}^b×(\boldsymbol{\omega}_{ib}^b×\mathbf{r}_N^b )) \end{matrix} \right] =\mathbf{J}\left[ \begin{matrix} \boldsymbol{\omega} ̇_{ib}^b \\ \mathbf{f}^b \end{matrix} \right]+\mathbf{M}\)     (5)

식 (5)를 각가속도와 비력에 대해 정리하면 식 (6)이 된다.

\(\left[ \begin{matrix} \boldsymbol{\omega} ̇_{ib}^b \\ \mathbf{f}^b \end{matrix} \right]=(\mathbf{J}^{\rm T} \mathbf{J})^{-1} \mathbf{J}^{\rm T} (\mathbf{Y}-\mathbf{M})=\left[ \begin{matrix} \mathbf{T}_{\boldsymbol{\omega} ̇} \\ \mathbf{T}_{\rm f} \end{matrix} \right] (\mathbf{Y}-\mathbf{M})\)                                             (6)

식 (6)에서 구한 비력 \(\mathbf{f}^b=\mathbf{T}_{\rm f} (\mathbf{Y}-\mathbf{M})\)를 식 (2)의 속도 미분방정식에 대입하면, 식 (7)과 같은 GF-INS의 속도 미분방정식을 얻을 수 있다.

\(\mathbf{v} ̇^n=\mathbf{C}_b^n [\mathbf{T}_{\rm f} (\mathbf{Y}-\mathbf{M})]-(2\boldsymbol{\Omega}_{ie}^n+\boldsymbol{\Omega}_{en}^n ) \mathbf{v}^n+\mathbf{g}_l^n\)                                                     (7)

따라서, GF-INS의 항법방정식은 식 (6)에서 각속도 미분방정식 \(\boldsymbol{\omega} ̇_{ib}^b=\mathbf{T}_{\boldsymbol{\omega} ̇} (\mathbf{Y}-\mathbf{M})\), 식 (1)의 자세 미분방정식, 식 (7)의 속도미분방정식, 그리고, 식 (3)의 위치 미분방정식으로 이루어진다.

2.2. GF-INS 오차모델

GF-INS 오차모델은 식 (6), (1), (7), (3)의 GF-INS 항법 방정식으로부터 구할 수 있다 (Klein 2015, Nusbaum & Klein 2017). 오차모델은 psi 각 오차모델과 phi 각 오차모델로 구분하는데, 본 논문에서는 실시간 항법 수행에 쓰이는 psi 각 오차모델을 이용한다. 식 (6)에 나타낸 각가속도와 비력을 섭동 (perturbation) 시키면, 식 (8)이 된다.

\(\left[ \begin{matrix} \boldsymbol{\omega} ̇_{ib}^b+\delta \boldsymbol{\omega} ̇_{ib}^b \\ \mathbf{f}^b+\delta \mathbf{f}^b \end{matrix} \right]=\left[ \begin{matrix} \mathbf{T}_{\dot\omega} \\ \mathbf{T}_{\rm f} \end{matrix} \right](\mathbf{Y}+\delta \mathbf{Y})-\left[ \begin{matrix} \mathbf{T}_{\dot\omega} \\ \mathbf{T}_{\rm f} \end{matrix} \right](\mathbf{M}+\delta \mathbf{M})\)                                   (8)

여기서 \(\mathbf{M}+\delta \mathbf{M}\)은 식 (9)와 같다.

\(\mathbf{M}+\delta \mathbf{M}=\left[ \begin{matrix} \mathbf{d}_1^{\rm T} (\boldsymbol{\Omega}_{ib}^b+\delta \boldsymbol{\Omega}_{ib}^b )(\boldsymbol{\Omega}_{ib}^b+\delta \boldsymbol{\Omega}_{ib}^b ) \mathbf{r}_1^b \\ ⋮ \\ \mathbf{d}_N^{\rm T} (\boldsymbol{\Omega}_{ib}^b+\delta \boldsymbol{\Omega}_{ib}^b )(\boldsymbol{\Omega}_{ib}^b+\delta \boldsymbol{\Omega}_{ib}^b ) \mathbf{r}_N^b \end{matrix} \right]\)                                                             (9)

식 (9)의 i번째 행을 다시 정리하면, 식 (10)이 된다.

\(\mathbf{d}_i^T (\boldsymbol{\Omega}_{ib}^b+\delta\boldsymbol{\Omega}_{ib}^b )(\boldsymbol{\Omega}_{ib}^b+\delta\boldsymbol{\Omega}_{ib}^b ) r_i^b \)                                                                                                                                                              
\(=\mathbf{d}_i^T \left[ \begin{matrix} (\omega_{ib,x}^b+\delta\omega_{ib,x}^b )(\omega_{ib,y}^b+\delta\omega_{ib,y}^b ) r_{i,y}^b+(\omega_{ib,x}^b+\delta\omega_{ib,x}^b )(\omega_{ib,z}^b+\delta\omega_{ib,z}^b ) r_{i,z}^b-((\omega_{ib,y}^b+\delta{\omega_{ib,y}^b}^2+(\omega_{ib,z}^b+\delta{\omega_{ib,z}^b}^2 ) r_{i,x}^b \\ (\omega_{ib,x}^b+\delta\omega_{ib,x}^b )(\omega_{ib,y}^b+\delta\omega_{ib,y}^b ) r_{i,x}^b+(\omega_{ib,y}^b+\delta\omega_{ib,y}^b )(\omega_{ib,z}^b+\delta\omega_{ib,z}^b ) r_{i,z}^b-((\omega_{ib,x}^b+\delta{\omega_{ib,x}^b }^2+(\omega_{ib,z}^b+\delta{\omega_{ib,z}^b}^2 ) r_{i,y}^b \\ (\omega_{ib,x}^b+\delta\omega_{ib,x}^b )(\omega_{ib,z}^b+\delta\omega_{ib,z}^b ) r_{i,x}^b+(\omega_{ib,y}^b+\delta\omega_{ib,y}^b )(\omega_{ib,z}^b+\delta\omega_{ib,z}^b ) r_{i,y}^b-((\omega_{ib,x}^b+\delta{\omega_{ib,x}^b }^2+(\omega_{ib,y}^b+\delta{\omega_{ib,y}^b }^2 ) r_{i,z}^b \end{matrix} \right]\)                              
\(=\mathbf{d}_i^T \left[ \begin{matrix} (\omega_{ib,x}^b+\delta\omega_{ib,x}^b )(\omega_{ib,y}^b+\delta\omega_{ib,y}^b ) r_{i,y}^b+(\omega_{ib,x}^b+\delta\omega_{ib,x}^b )(\omega_{ib,z}^b+\delta\omega_{ib,z}^b ) r_{i,z}^b-((\omega_{ib,y}^b+\delta{\omega_{ib,y}^b}^2+(\omega_{ib,z}^b+\delta{\omega_{ib,z}^b}^2 ) r_{i,x}^b \\ (\omega_{ib,x}^b+\delta\omega_{ib,x}^b )(\omega_{ib,y}^b+\delta\omega_{ib,y}^b ) r_{i,x}^b+(\omega_{ib,y}^b+\delta\omega_{ib,y}^b )(\omega_{ib,z}^b+\delta\omega_{ib,z}^b ) r_{i,z}^b-((\omega_{ib,x}^b+\delta{\omega_{ib,x}^b }^2+(\omega_{ib,z}^b+\delta{\omega_{ib,z}^b}^2 ) r_{i,y}^b \\ (\omega_{ib,x}^b+\delta\omega_{ib,x}^b )(\omega_{ib,z}^b+\delta\omega_{ib,z}^b ) r_{i,x}^b+(\omega_{ib,y}^b+\delta\omega_{ib,y}^b )(\omega_{ib,z}^b+\delta\omega_{ib,z}^b ) r_{i,y}^b-((\omega_{ib,x}^b+\delta{\omega_{ib,x}^b }^2+(\omega_{ib,y}^b+\delta{\omega_{ib,y}^b }^2 ) r_{i,z}^b \end{matrix} \right] \)                              
\(+\mathbf{d}_i^T \left[ \begin{matrix} \omega_{ib,y}^b r_{i,y}^b+\omega_{ib,z}^b r_{i,z}^b&\omega_{ib,x}^b r_{i,y}^b-2\omega_{ib,y}^b r_{i,x}^b&\omega_{ib,x}^b r_{i,z}^b-2\omega_{ib,z}^b r_{i,x}^b \\ \omega_{ib,y}^b r_{i,x}^b-2\omega_{ib,x}^b r_{i,y}^b&\omega_{ib,x}^b r_{i,x}^b+\omega_{ib,z}^b r_{i,z}^b&\omega_{ib,y}^b r_{i,z}^b-2\omega_{ib,z}^b r_{i,y}^b \\ \omega_{ib,z}^b r_{i,x}^b-2\omega_{ib,x}^b r_{i,z}^b&\omega_{ib,z}^b r_{i,y}^b-2\omega_{ib,y}^b r_{i,z}^b&\omega_{ib,x}^b r_{i,x}^b+\omega_{ib,y}^b r_{i,y}^b \end{matrix} \right] \left[ \begin{matrix} \delta\omega_{ib,x}^b \\ \delta\omega_{ib,y}^b \\ \delta\omega_{ib,z}^b \end{matrix} \right]\)                   (10)

여기서 \(\omega_{ib,x}^b,\omega_{ib,y}^b,\omega_{ib,z}^b\)는 각각 x, y, z축의 각속도, \(\delta\omega_{ib,x}^b,\delta\omega_{ib,y}^b,\delta\omega_{ib,z}^b\)는 각각 x, y, z축의 각속도 오차, \(r_{i,x}^b,r_{i,y}^b,r_{i,z}^b\)는 각각 위치의 x, y, z축 성분을 의미한다. 식 (10)에서 각속도 오차에 대한 항만 남기면, \(\delta \mathbf{M}=\left[ \begin{matrix} \delta m_1 & ⋯ & \delta m_N \end{matrix} \right]^T\)의 i번째 항 \(\delta m_1\)는 식 (11)이 된다.

\(\delta m_i=\left[ \begin{matrix} d_{i,x}&d_{i,y}&d_{i,z} \end{matrix} \right] \left[ \begin{matrix} \omega_{ib,y}^b r_{i,y}^b+\omega_{ib,z}^b r_{i,z}^b&\omega_{ib,x}^b r_{i,y}^b-2\omega_{ib,y}^b r_{i,x}^b&\omega_{ib,x}^b r_{i,z}^b-2\omega_{ib,z}^b r_{i,x}^b \\ \omega_{ib,y}^b r_{i,x}^b-2\omega_{ib,x}^b r_{i,y}^b&\omega_{ib,x}^b r_{i,x}^b+\omega_{ib,z}^b r_{i,z}^b&\omega_{ib,y}^b r_{i,z}^b-2\omega_{ib,z}^b r_{i,y}^b \\ \omega_{ib,z}^b r_{i,x}^b-2\omega_{ib,x}^b r_{i,z}^b&\omega_{ib,z}^b r_{i,y}^b-2\omega_{ib,y}^b r_{i,z}^b&\omega_{ib,x}^b r_{i,x}^b+\omega_{ib,y}^b r_{i,y}^b \end{matrix} \right] \delta \boldsymbol{\omega}_{ib}^b\)                
\(=\left[ \begin{matrix} a_{i,1}&a_{i,2}&a_{i,3} \end{matrix} \right] \delta \boldsymbol{\omega}_{ib}^b\)                                                                                                                                                   (11)

따라서, 식 (8)에 나타낸 \(\delta \mathbf{M}\)은 식 (12)가 된다.

\(\delta \mathbf{M}=\left[ \begin{matrix} a_{1,1}&a_{1,2}&a_{1,3}\\⋮&⋮&⋮\\a_{N,1}&a_{N,2}&a_{N,3} \end{matrix} \right]\delta \boldsymbol{\omega}_{ib}^b=\mathbf{A}\delta \boldsymbol{\omega}_{ib}^b\)                                                           (12)

위의 결과를 이용하여 식 (8)의 각가속도 오차와 비력 오차를 구하면, 각각 식 (13)과 (14)가 된다.

\(\delta \boldsymbol{\omega ̇}_{ib}^b=\mathbf{T}_\boldsymbol{\omega} ̇ (\delta \mathbf{Y}-\mathbf{A}\delta \boldsymbol{\omega}_{ib}^b )\)                                                                                         (13)

\(\delta \mathbf{f}^b=\mathbf{T}_{\rm f} (\delta \mathbf{Y}-\mathbf{A}\delta \boldsymbol{\omega}_{ib}^b )\)                                                                                          (14)

마찬가지로 위치, 속도, 자세 오차 모델을 구하면, 식 (15-17)이 된다.

\(\delta \boldsymbol{\Psi ̇ }=-\mathbf{C}_b^n \delta \boldsymbol{\omega}_{ib}^b-\boldsymbol{\Omega}_{in}^n \delta \boldsymbol{\Psi}\)                                                                                      (15)

\(\delta \mathbf{v ̇}^n=\mathbf{C}_b^n (\mathbf{T}_{\mathbf{f}} \delta \mathbf{Y}-\mathbf{T}_{\mathbf{f}} \mathbf{A}\delta \boldsymbol{\omega}_{ib}^b )+[\mathbf{f}^n×]\delta \boldsymbol{\Psi}-(2\boldsymbol{\Omega}_{ie}^n+\boldsymbol{\Omega}_{en}^n )\delta \mathbf{v}^n\)                                        (16)

\(\delta \mathbf{p} ̇^n=-\boldsymbol{\Omega}_{en}^n \delta \mathbf{p}^n+\delta \mathbf{v}^n\)                                                                                          (17)

여기서 \(\delta \boldsymbol{\Psi }\)는 자세오차, \(\delta \boldsymbol{\omega}_{ib}^b\)는 각속도 오차, \(\delta \mathbf{v}^n\)는 속도 오차, \(\delta \mathbf{p}^n\)는 위치 오차를 나타낸다.

전술한 GF-INS의 오차모델과 INS의 오차모델을 비교하면, INS의 각속도 오차는 자이로 오차이며, GF-INS의 각속도 오차는 각가속도 오차가 누적된 값이다. 따라서, GF-INS의 각속도 오차는 시간이 지남에 따라 INS의 각속도 오차보다 한 차수가 높은 속도로 발산하며, 그 결과 자세 오차도 더 큰 것을 예상할 수 있다.

2.3. 통합 칼만필터

통합 칼만필터를 구성하기 위한 상태방정식은 식 (18)과 같다.

\(\delta \mathbf{x} ̇= \mathbf{F} \delta \mathbf{x}+\mathbf{w}\)(18)

여기서 \(\delta \mathbf{x}\)는 상태 변수, \(\mathbf{F}\)는 시스템 행렬, \(\mathbf{w}\)는 공분산 행렬이 \(\mathbf{Q}\)인 공정 잡음 벡터를 의미한다. 칼만필터 구성을 위해서 GF-INS 오차모델인 식 (13), (15), (16), (17)을 이용한다. 가속도계 오차 \((\delta \mathbf{Y})\)는 식 (19)와 같이 모델링한다.

\(\delta \mathbf{Y} ̇= \mathbf{b_A} + \mathbf{w_A}\)(19)

여기서 \(\mathbf{b_A}\)는 가속도계 바이어스, \(\mathbf{w_A}\)는 가속도계 잡음을 의미한다. 식 (18)에서 상태변수는 식 (20)과 같다.

\(\delta \mathbf{x}=\left[ \begin{matrix} \delta \mathbf{p}^{nT} & \delta \mathbf{v}^{nT} & \delta \boldsymbol{\Psi}^T & \delta \boldsymbol{\omega}_{ib}^{bT} & \mathbf{b_A}^T \end{matrix} \right]^T\)(20)

여기서 위치 오차 \(\delta \mathbf{p}^{n}\), 속도 오차 \(\delta \mathbf{v}^{n}\), 자세 오차 \(\delta \mathbf{\Psi}\), 각속도 오차 \(\delta \boldsymbol{\omega}_{ib}^{b}\), 바이어스 오차 \(\mathbf{b_A}\)는 각각 식 (21-25)와 같다.

\(\delta \mathbf{p}^n=\left[ \begin{matrix} \delta p_N^n&\delta p_E^n&\delta p_D^n \end{matrix} \right]^T\)                                                                             (21)

\(\delta \mathbf{v}^n=\left[ \begin{matrix} \delta v_N^n&\delta v_E^n&\delta v_D^n \end{matrix} \right]^T\)                                                                              (22)

\(\delta \boldsymbol{\Psi}=\left[ \begin{matrix} \delta \alpha&\delta \beta&\delta \gamma \end{matrix} \right]^T\)                                                                                       (23)

\(\delta \boldsymbol{\omega}_{ib}^b=\left[ \begin{matrix} \delta \omega_{ib,x}^b&\delta \omega_{ib,y}^b&\delta \omega_{ib,z}^b \end{matrix} \right]^T\)                                                                   (24)

\(\mathbf{b_A} =\left[ \begin{matrix} b_{A1}&⋯&b_{AN} \end{matrix} \right]^T\)                                                                                  (25)

여기서 \(\delta p_N^n\)는 north축 위치 오차, \(\delta p_E^n\)는 east축 위치 오차, \(\delta p_D^n\)는 down축 위치 오차, \(\delta \alpha\)는 roll 오차, \(\delta \beta\)는 pitch 오차, \(\delta \gamma\)는 yaw 오차, \(\delta \omega_{ib,x}^b\)는 x축의 각속도 오차, \(\delta \omega_{ib,y}^b\)는 y축의 각속도 오차, \(\delta \omega_{ib,z}^b\)는 z축의 각속도 오차, \(b_{A1}\)는 i번째 가속도계 바이어스를 의미한다. 시스템 행렬 \(\mathbf{F}\)는 식 (26)과 같다.

\(\mathbf{F}=\left[ \begin{matrix} -\boldsymbol{\Omega}_{en}^n&\mathbf{I}_{3×3}&0_{3×3}&0_{3×3}&0_{3×N} \\ 0_{3×3}&-(2\boldsymbol{\Omega}_{ie}^n+\boldsymbol{\Omega}_{en}^n )&[\mathbf{f}^n×]&-\mathbf{C}_b^n \mathbf{T_f} \mathbf{A}&\mathbf{C}_b^n \mathbf{T_f} \\ 0_{3×3}&0_{3×3}&-\boldsymbol{\Omega}_{in}^n&-\mathbf{C}_b^n&0_{3×N} \\ 0_{3×3}&0_{3×3}&0_{3×3}&-\mathbf{T}_\omega ̇ \mathbf{A}&\mathbf{T}_\omega ̇ \\ 0_{N×3}&0_{N×3}&0_{N×3}&0_{N×3}&0_{N×N} \end{matrix} \right]\)                                               (26)

GPS 수신기의 위치와 속도만을 이용할 경우의 측정방정식은 식 (27)과 같다.

\(\delta \mathbf{z}_1=\left[ \begin{matrix} \mathbf{p}_{GF-INS}^n-\mathbf{p}_{GPS}^n \\ \mathbf{v}_{GF-INS}^n-\mathbf{v}_{GPS}^n \end{matrix} \right] +\boldsymbol{\upsilon}_1=\left[ \begin{matrix} \mathbf{I}_{6×6}&0_{6×6}&0_{6×N} \end{matrix} \right] \delta \mathbf{x}+\boldsymbol{\upsilon}_1\)                                               (27)

여기서 \(\mathbf{p}_{GF-INS}^n\)는 GF-INS 위치, \(\mathbf{p}_{GPS}^n\)는 GPS 수신기 위치, \(\mathbf{v}_{GF-INS}^n\)는 GF-INS 속도, \(\mathbf{v}_{GPS}^n\)는 GPS 수신기 속도, \(\boldsymbol{\upsilon}_1\)은 첫번째 측정치의 잡음을 의미한다. GPS 수신기의 위치, 속도와 자세 차이가 입력될 경우의 측정방정식은 식 (28)과 같다.

\(\delta \mathbf{z}_2=\left[ \begin{matrix} \mathbf{p}_{GF-INS}^n-\mathbf{p}_{GPS}^n \\ \mathbf{v}_{GF-INS}^n-\mathbf{v}_{GPS}^n \\ \boldsymbol{\Psi}_{GF-INS}-\boldsymbol{\Psi}_{GPS} \end{matrix} \right]+\boldsymbol{\upsilon}_2 \)                                                                                                                     
\( = \left[ \begin{matrix} \mathbf{I}_{6×6}&0_{3×1}&0_{3×1}&0_{3×1}&0_{3×3}&0_{6×N} \\ 0_{1×6}&-{\cos ⁡ \gamma \over \cos ⁡\beta} &- {\sin ⁡\gamma \over \cos ⁡\beta} &0&0_{1×3}&0_{1×N} \\ 0_{1×6}&\sin ⁡\gamma & -\cos ⁡\gamma&0&0_{1×3}&0_{1×N} \\ 0_{1×6}&-\tan ⁡\beta \cos ⁡\gamma&-\tan ⁡\beta \sin ⁡\gamma&-1&0_{1×3}&0_{1×N} \end{matrix} \right] \delta \mathbf{x} +\boldsymbol{\upsilon}_2\)                         (28)

여기서 \({\Psi}_{GF-INS}\)는 GF-INS의 오일러 각, \(\boldsymbol{\Psi}_{GPS}\)는 GPS 수신기에서 출력된 오일러 각, \(\boldsymbol{\upsilon}_2\)은 두번째 측정치의 잡음을 의미한다.

단일 안테나를 사용하는 일반적인 GPS 수신기는 항체의 위치와 속도를 출력한다 (Titterton & Weston 2004). 그리고, 다중 안테나를 사용하는 GPS 수신기는 항체의 위치, 속도와 자세를 제공한다. 이러한 GPS 수신기는 Septentrio (2009)사의 AsteRx2eH 수신기와 Trimble사의 TANS Vector GPS가 있다 (Li et al. 2004). 관련 연구 결과로 Park (1996)과 Moon (2014)이 다중 안테나를 사용한 GPS 수신기로 자세를 추정한 것을 볼 수 있다.

3. 제안한 알고리즘의 성능평가

제안한 GPS/GF-INS 통합항법 시스템의 성능을 보기 위하여 각속도가 일정한 circling 궤적과 각속도가 변하는 spiral 궤적에 대한 항법 결과를 확인하였다. 출력율을 높인 GPS 수신기의 위치, 속도, 자세가 보조 신호인 제안한 GPS/GF-INS 통합항법 시스템의 성능은 시뮬레이션을 통해 확인하였다. Fig. 4와 같은 가속도계 6개로 구성된 Cube 구조의 GF-IMU를 대상으로 시뮬레이션을 수행하였는데, 가속도계 6개의 위치는 원점으로부터 만큼 떨어진 세 축 위에 각각 대칭되도록 두었으며, 각 면에 45도 대각이 측정방향이다. 이 경우 6개의 가속도계 바이어스는 상태변수 식 (20)의 \(\mathbf{b_A}\)이며, \(\mathbf{T_\dot{\boldsymbol{\omega}}}\)와 \(\mathbf{T_f}\)는 각각 식 (29)와 (30)이 된다.

Fig. 4. GF-IMU for performance evaluation.

\(\mathbf{T_\dot{\boldsymbol{\omega}}}={1 \over 2\sqrt{2} l} \left[ \begin{matrix} 1&-1&0&0&1&-1 \\ -1&0&1&-1&0&-1 \\ 0&1&-1&-1&1&0 \end{matrix} \right]\)                                                             (29)

\(\mathbf{T_f}={1 \over 2\sqrt{2}} \left[ \begin{matrix} 1&1&0&0&-1&-1 \\ 1&0&1&-1&0&1 \\ 0&1&1&1&1&0 \end{matrix} \right]\)                                                                     (30)

\(\mathbf{T_\dot{\boldsymbol{\omega}}} \mathbf{A}\)와 \(\mathbf{T_f} \mathbf{A}\)는 각각 식 (31), (32)와 같다.

\(T_{\dot{\omega}} A=0_{3×3}\)                                                                                                               (31)

\(\mathbf{T_f A}=-l \left[ \begin{matrix} 0&ω_{ib,z}^b&ω_{ib,y}^b \\ ω_{ib,z}^b&0&ω_{ib,x}^b \\ ω_{ib,y}^b&ω_{ib,x}^b&0 \end{matrix} \right]=-l\mathbf{W}\)                                                                 (32)

따라서, 시스템 행렬 \(\mathbf{F}\)는 식 (26)에서 N이 6인 식 (33)이 된다.

\(\mathbf{F}=\left[ \begin{matrix} -\boldsymbol{\Omega}_{en}^n&\mathbf{I}_{3×3}&0_{3×3}&0_{3×3}&0_{3×6} \\ 0_{3×3} & -(2\boldsymbol{\Omega}_{ie}^n+\boldsymbol{\Omega}_{en}^n) & [\mathbf{f}^n×]&l\mathbf{C}_b^n \mathbf{W}&\mathbf{C}_b^n \mathbf{T_f} \\ 0_{3×3}&0_{3×3}&-\boldsymbol{\Omega}_{in}^n&-\mathbf{C}_b^n&0_{3×6} \\ 0_{3×3}&0_{3×3}&0_{3×3}&0_{3×3}&\mathbf{T}_{\dot\omega} \\ 0_{6×3}&0_{6×3}&0_{6×3}&0_{6×3}&0_{6×6} \end{matrix} \right]\)                                                 (33)

GPS 수신기의 위치, 속도만을 이용할 경우, 측정방정식은 식 (27)에서 N이 6인 것이고, GPS 수신기의 위치, 속도, 자세를 이용할 경우에는 식 (28)에서 N이 6인 것이 된다.

원점으로부터 가속도계 위치까지의 거리 l은 0.1 m, 가속도계는 전술급 (tactical) IMU에서 사용하는 Honeywell사 제품인 RBA-500을 기준으로 두었다. Table 1에는 RBA-500의 사양, Table 2에는 GPS 수신기의 사양을 나타내었다.

Table 1. RBA500 accelerometer specification.

Velocity random walk (1 σ)	Bias (1 σ)	Output rate
0.0066 m/s/\( \sqrt{\rm hr} \)	1 mg	100 Hz

 

Table 2. GPS receiver specification.

Pseudorange error (1 σ)	Pseudorange rate error (1 σ)	Attitude error (1 σ)	Clock bias (1 σ)	Clock drift (1 σ)	Output rate
1.5 m	0.05 m/s	0.5 deg	\(3.1623 \times 10^{-10} {\rm m}\)	\(6.2832 \times 10^{-10} {\rm m}\)	1 Hz, 10 Hz

 

3.1 Circling 궤적에 대한 성능평가

Fig. 5a에는 항체의 circling 궤적, Fig. 5b에는 시간에 따른 항체의 위치, 속도, 자세를 나타내었다. Fig. 5b를 보면, Roll과 Pitch는 일정하며, Yaw만 변하며, Yaw의 속도가 일정한 것을 알 수 있다.

Fig. 5. (a) Trajectory of the vehicle, (b) Position, velocity, and attitude of the vehicle.

Fig. 6. Simulation result with 1 Hz GPS receiver for circling trajectory. (a) Position error, (b) Velocity error, (c) Attitude error, (d) Angular velocity error

Fig. 7. Simulation result with 10 Hz GPS receiver for circling trajectory. (a) Position error, (b) Velocity error, (c) Attitude error, (d) Angular velocity error

Fig. 8. Simulation result with multi-antenna 1 Hz GPS receiver for circling trajectory. (a) Position error, (b) Velocity error, (c) Attitude error, (d) Angular velocity error

Fig. 9. Simulation result with multi-antenna 10 Hz GPS receiver for circling trajectory. (a) Position error, (b) Velocity error, (c) Attitude error, (d) Angular velocity error

GPS 수신기가 위치와 속도를 제공하고, 출력율이 1 Hz인 경우, GPS/GF-INS 통합항법 시스템의 위치, 속도, 자세, 각속도 오차를 Fig. 6에 나타내었으며, 출력율이 10 Hz인 경우의 위치, 속도, 자세, 각속도 오차를 Fig. 7에 나타내었다. GPS 수신기가 위치, 속도, 자세를 제공하고, 출력율이 1 Hz인 경우, GPS/GF-INS 통합항법 시스템의 위치, 속도, 자세, 각속도 오차를 Fig. 8, 출력율이 10 Hz인 경우의 위치, 속도, 자세, 각속도 오차를 Fig. 9에 나타내었다. Table 3에는 Figs. 6-9에 나타낸 항법 결과에 대하여 위치, 속도, 자세, 각속도의 Root Mean Square Error (RMSE)를 나타내었다. Table 3을 보면, GPS 수신기가 위치와 속도만 제공할 경우에는 출력율을 높이면, 위치와 속도 추정 성능이 매우 좋아지지만, 자세 오차와 각속도 오차는 매우 큰 것을 알 수 있다. 반면에, GPS 수신기에서 위치, 속도와 자세를 제공할 경우에는 출력율을 높이면, 자세 오차와 각속도 오차가 크게 줄어드는 것을 알 수 있다.

Table 3. RMSE of navigation results.

GPS receiver outputs		Position, velocity		Position, velocity, attitude
Output rate (Hz)		1	10	1	10
Position (m)	North  East  Down	0.50 0.49 1.73	0.10 0.08 0.17	0.22 0.29 0.29	0.06 0.08 0.17
Velocity (m/s)	North  East  Down	0.40 0.27 1.05	0.01 0.02 0.02	0.07 0.08 0.14	0.01 0.01 0.02
Attitude (deg)	Roll  Pitch  Yaw	107.80 34.84 150.86	4.80 24.75 14.85	0.91 2.40 2.07	0.09 0.14 0.13
Angular velocity (deg/s)	x  y  z	4.32 11.04 10.38	3.05 0.83 1.83	0.72 0.30 2.12	0.19 0.09 0.27

 

3.2 Spiral 궤적에 대한 성능평가

Fig. 10a에는 항체의 spiral 궤적, Fig. 10b에는 시간에 따른 항체의 위치, 속도, 자세를 나타내었다. Fig. 10b에서 항체의 자세를 보면, Pitch와 Yaw가 변하며, Yaw의 속도가 일정한 것을 알 수 있다.

Fig. 10. (a) Trajectory of the vehicle, (b) Position, velocity, and attitude of the vehicle.

GPS 수신기가 위치와 속도를 제공하고, 출력율이 1 Hz인 경우의 위치, 속도, 자세, 각속도 오차를 Fig. 11에 나타내었으며, 출력율이 10 Hz인 경우의 위치, 속도, 자세, 각속도 오차를 Fig. 12에 나타내었다. GPS 수신기가 위치, 속도, 자세를 제공하고, 출력율이 1 Hz인 경우의 위치, 속도, 자세, 각속도 오차를 Fig. 13에 나타내었으며, 출력율이 10 Hz인 경우의 위치, 속도, 자세, 각속도 오차를 Fig. 14에 나타내었다. Table 4에는 Figs. 11-14에 나타낸 항법 결과에 대하여 위치, 속도, 자세, 각속도의 RMSE를 나타내었다. Table 4을 보면, 제안한 통합항법 시스템처럼 자세 측정치를 추가하고, GPS 수신기의 출력율을 높이면, 자세 오차와 각속도 오차가 줄어드는 것을 알 수 있다. 그리고, GPS 수신기가 1 Hz일 때, 자세 측정치를 추가한 결과를 보면, Table 3의 결과에 비해 위치, 속도, 자세, 각속도 오차가 크게 줄어든 것을 알 수 있다.

Fig. 11. Simulation result with 1 Hz GPS receiver for spiral trajectory. (a) Position error, (b) Velocity error, (c) Attitude error, (d) Angular velocity error

Fig. 12. Simulation result with 10 Hz GPS receiver for spiral trajectory. (a) Position error, (b) Velocity error, (c) Attitude error, (d) Angular velocity error

Fig. 13. Simulation result with multi-antenna 1 Hz GPS receiver for spiral trajectory. (a) Position error, (b) Velocity error, (c) Attitude error, (d) Angular velocity error

Fig. 14. Simulation result with multi-antenna 10 Hz GPS receiver for spiral trajectory. (a) Position error, (b) Velocity error, (c) Attitude error, (d) Angular velocity error

Table 4. RMSE of navigation results.

GPS receiver outputs		Position, velocity		Position, velocity, attitude
Output rate (Hz)		1	10	1	10
Position (m)	North  East  Down	26699.46 7920.32 80481.53	0.12 0.18 0.20	28580.52 5926.69 80645.61	0.09 0.20 0.22
Velocity (m/s)	North  East  Down	92.61 258.90 1048.72	0.02 0.02 0.01	462.74 195.52 1051.34	0.02 0.02 0.01
Attitude (deg)	Roll  Pitch  Yaw	80.73 43.64 102.85	17.28 22.78 65.37	110.09 42.52 103.37	0.11 0.14 0.18
Angular velocity (deg/s)	x  y  z	81.81 371.13 35.53	0.74 0.55 2.89	137.81 276.04 96.53	0.11 0.18 0.12

 

Figs. 6, 8, 11을 보면, 1 Hz GPS 수신기를 사용할 경우, GPS/GF-INS 통합항법 결과가 발산하는 것을 알 수 있다. 식 (17)에 나타낸 GF-INS의 위치 오차 \((\delta \mathbf{p}^n)\)는 t⁴에 비례하는 반면에, INS의 위치 오차는 t3에 비례한다. 따라서, INS보다 GF-INS의 항법 오차가 더 빠르게 증가하므로, 이를 막기 위하여 출력율을 높이는 것이 필요함을 알 수 있다. 그리고, 정지 시의 GPS/GF-INS 통합 칼만필터의 가관측성을 확인하면, GPS 수신기의 출력이 위치와 속도인 경우에는 가관측하지 않지만, 자세까지 추가하면, 가관측한 것을 쉽게 확인할 수 있다. 그런데, 본 논문에서 대상으로 하는 궤적은 정지한 것이 아니므로, GPS 출력이 위치와 속도인 경우에는 발산하는 것을 예상할 수 있다. 이러한 특성은 정지 시 정렬을 하기 위하여 GPS/GF-INS 통합 칼만필터를 구성한 Nusbaum & Klein (2017)의 가관측성 분석 결과에서도 볼 수 있다. 즉, 빠르게 증가하는 GF-INS의 항법 오차 보정을 위하여 빠른 출력율의 GPS 수신기가 필요하며, 위치와 속도, 자세 측정치가 모두 필요함을 알 수 있다.

일반적으로 GPS 수신기의 출력율은 1 Hz로 알려져 있으나, Table 2에 나타낸 GPS 수신기의 출력율은 10 Hz이다. 최근 출시된 u-blox사의 NEO-7, NEO-M8 수신기의 출력율인 10 Hz, Swift Navigation사의 Piksi Multi 수신기의 출력율도 10 Hz이므로, 이것은 쉽게 이용가능한 것임을 알 수 있다 (u-blox 2015, 2021, Swift Navigation 2019).

4. CONCLUSIONS

본 논문에서는 위치, 속도뿐만 아니라 자세 정보까지 제공하는 GPS와 GPS의 출력율을 높이는 방법을 이용하는 GPS/GF-INS 통합항법 알고리즘을 제안하고, 성능을 평가한 결과를 제시하였다. GF-INS의 항법 방정식과 오차모델을 정리하고, 정육면체의 각 면에 가속도계를 배치한 구조에 대하여 항체가 circling 궤적과 spiral 궤적으로 운동할 때, 제안한 통합항법의 유효성을 확인하였다. 성능평가 결과 GPS 수신기의 위치와 속도만을 보조 신호로 사용할 경우에 비하여 GPS 수신기의 위치, 속도, 자세를 보조 신호로 사용한 제안한 통합항법 시스템이 더 정확한 항법 결과를 얻을 수 있다. 그리고, 보조 신호인 GPS의 출력율을 높일 경우, 더 정확한 항법 결과를 얻을 수 있다. 추후에는 제안한 항법을 이용하는 GPS/GF-INS 통합항법 시스템을 항체에 장착하여 실험을 수행하고, 여러가지 궤적에 대하여 가관측성을 분석할 예정이다.