Bearing-only Localization of GNSS Interference using Iterated Consider Extended Kalman Filter

Abstract

In this paper, the Iterated Consider Extended Kalman Filter (ICEKF) is proposed for bearing-only localization of GNSS interference to improve the estimation performance and filter consistency. The ICEKF is an extended version of Consider KF (CKF) for Iterated EKF (IEKF) to consider an effect of bearing measurement bias error to filter covariance. The ICEKF can mitigate the EKF divergence problem which can occur when linearizing the nonlinear bearing measurement by a large initial state error. Also, it can mitigate filter inconsistency problem of EKF and IEKF which can occur when a weakly observable bearing measurement bias error state is not included in filter state vector. The simulation result shows that the localization error of the ICEKF is smaller than the EKF and IEKF, and the Root Mean Square (RMS) estimation error of ICEKF matches the covariance of filter.

1. INTRODUCTION

Global Positioning System (GPS)로 대표되는 전지구 위성항법시스템(Global Navigation Satellite System, GNSS)은 의도적 또는 비의도적 전파교란에 취약한 특성이 있기 때문에 GNSS 수신기를 운용하는 사용자에게 GNSS 교란에 대한 감시 정보를 제공하는 GNSS 교란감시 시스템의 역할이 중요하다. GNSS 교란감시 시스템은 교란탐지, 신호특성분석, 교란원 위치탐지 및 교e란영향 분석 등 사용자에게 필요한 다양한 정보를 제공할 수 있는데, 교란에 의한 영향을 분석하기 위해서는 교란원의 위치탐지가 선행되어야 한다. 교란원의 위치탐지는 다수의 센서에서 측정된 교란신호의 도래각(Angle Of Arrival, AOA), 도래 시간 차이(Time Difference Of Arrival, TDOA), 도래 주파수 차이(Frequency Difference Of Arrival, FDOA), 및 도래 전력 차이(Power Difference Of Arrival, PDOA) 등의 측정치를 이용하여 이루어질 수 있다 (Hoye 2010). 본 논문에서는 배열안테나를 이용하여 교란신호의 방위각 도래각을 측정하는 센서들로 구성되는 지상기반 GNSS교란감시 시스템에서 칼만필터를 이용한 위치탐지 기법을 다룬다.

방위각 도래각의 측정식은 센서의 위치와 교란원의 위치에 대한 비선형함수로 표현되기 때문에 칼만필터로 위치를 추정하는 경우 기본적으로 확장칼만필터(Extended Kalman Filter, EKF)로 구현하며, EKF는 상태변수와 오차공분산의 초기값으로부터 연속된 측정치에 대한 측정치 갱신을 통해 위치를 추정하는데, 상태변수의 초기값이 실제 교란원의 위치와 큰 차이가 나는 경우 비선형 측정식에 대한 선형화오차에 의해 위치 추정오차가 크게 발생하거나 필터가 발산할 수 있다 (Tully et al. 2008). 특히, 도래각 기반의 위치추정 정확도는 측정치의 정확도와 교란원 및 센서들간의 배치형상에 따른 Dilution Of Precision (DOP)에 따라 결정되기 때문에 (Dempster 2006), DOP이 큰 위치에 존재하는 교란원인 경우 EKF가 발산할 수 있다. 이와 같이 비선형 측정식에서 상태변수의 초기값 오차가 큰 경우 발생하는 EKF의 발산을 억제하기 위해 동일한 측정치에 대해 측정치 갱신을 반복하여 선형화 오차에 의한 영향을 감소시키는 반복확장필터(Iterated Extended Kalman Filter, IEKF)를 적용할 수 있다(Havlik & Straka 2015).

칼만필터는 대상 시스템의 특성과 오차 요소를 모두 포함하여 모델링하는 경우 최적필터로 설계가 가능하다. 그러나 상태변수의 크기가 매우 커져 구현하기에 부적절한 경우, 또는 가관측성이 낮아 추정이 거의 되지 않는 상태변수나 항법오차 추정성능 향상에 미치는 영향이 미미한 상태변수가 있는 경우에는 차수가 감소된 준최적 칼만필터로 설계한다. 준최적 칼만필터를 구성하는 방법으로 불필요 상태변수를 단순히 제거하는 차수 감소필터(pure reduced order filter)와 불필요 상태변수를 추정하지는 않지만 오차특성을 공분산 계산에 반영하는 반영칼만필터(Consider Kalman Filter, CKF)가 있다 (Brown & Hwang 1997). 센서의 도래각 측정치에는 잡음오차 뿐만 아니라 바이어스형태의 오차가 포함될 수 있는데, 바이어스오차는 가관측하지만 거의 추정이 되지 않기 때문에 필터 상태변수에서 제외하고 칼만필터의 측정잡음에 특성을 반영하는 차수 감소 EKF 및 IEKF로 구현할 수 있지만, 단순 차수 감소필터는 상태변수 추정치와 칼만필터 공분산이 일치하지 않는 필터의 일관성(consistency)이 깨지는 문제가 발생할 수 있다.

본 논문에서는 교란신호의 방위각 도래각을 이용한 위치추정 기법으로 DOP이 나쁜 경우에 큰 초기 위치오차에 의한 EKF의 발산을 억제하고 도래각 측정치에 잡음오차 뿐만 아니라 바이어스 형태의 오차가 포함되는 경우에 바이어스 오차의 특성을 필터 공분산 계산에 반영할 수 있는 Iterated Consider Extended Kalman Filter (ICEKF)를 설계하고 시뮬레이션을 통해 ICEKF와 차수 감소 EKF 및 IEKF의 성능을 위치추정 정확도와 필터 일관성 관점에서 분석하였다.

2. ITERATED EXTENDED kALMAN FILTER WITH CONSIDER STATE

2.1 Extended Kalman Filter

본 논문에서는 정지된 전파교란원에 대해 방위각 도래각을 이용한 위치추정 문제를 다루고 있기 때문에 Eqs. (1-2)와 같이 이산시간에 대한 선형 시스템 모델과 비선형 측정모델을 정의한다.

\({\bf x}_{k+1}={\bf \Phi}_k{\bf x}_k+{\bf w}_k \)       (1)

\({\bf z}_k={\bf h}( {\bf x}_k )+{\bf v}_k\)       (2)

여기서, \({\bf\it x}_k\)는  \(n\times1 \)는 상태변수 벡터, \({\bf \Phi}_k\)는 \( n\times n\) 상태천이행렬, \({\bf z}_k\)는 \(m\times1\) 측정치 벡터, \({\bf\it h}({\bf\it x}_k)\)는 상태변수에 대한 비선형 측정식이며, 시스템 잡음 \({\bf w}_k\)와 측정잡음 \({\bf v}_k\)는 각각 공분산이 \(\mathbf{Q}\), \({\bf \it R}\)인 영평균 가우시안 잡음 시퀀스로 Eq. (3)을 만족한다.

\(E\left[{\bf w}_i{\bf w}_j^T\right]={\bf Q}\delta_{ij}\ ,\ \ E\left[{\bf v}_i{\bf v}_j^T\right]={\bf R}\delta_{ij}\ ,\ \ E\left[{\bf w}_i{\bf v}_j^T\right]=\mathbf{0} \)       (3)

여기서, \(\delta_{ij}\)는 Kronecker delta 함수이다.

Eq. (1)의 선형 시스템모델에 대한 칼만필터의 시간전파는 Eqs. (4-5)와 같다 (Simon 2006).

\({\hat{{\bf x}}}_k^-={\bf \Phi}_{k-1}{\hat{{\bf x}}}_{k-1}^+\)       (4)

\({\bf P}_k^-={\bf \Phi}_{k-1}{\bf P}_{k-1}^+{\bf \Phi}_{k-1}^T+{\bf Q}\)       (5)

여기서, \({\hat{{\bf x}}}_k\)는 상태변수 추정치, \({\bf P}_k\)는 상태변수 추정치에 대한 오차 공분산, 그리고 위첨자 -와 +는 각각 시간전파 및 측정치 갱신된 변수를 의미한다.

Eq. (2)의 비선형 측정모델에 대한 EKF의 측정치 갱신은 Eqs. (6-8)과 같다 (Simon 2006).

\({\bf K}_k={\bf P}_k^-{\bf H}_k^T\left[{\bf H}_k{\bf P}_k^-{\bf H}_k^T+{\bf R}\right]^{-1}\)       (6)

\({\hat{{\bf x}}}_k^+={\hat{{\bf x}}}_k^-+{\bf K}_k\left({\bf z}_k-{\bf h}\left({\hat{{\bf x}}}_k^-\right)\right)\)       (7)

\({\bf P}_k^+=\left[{\bf I}-{\bf K}_k{\bf H}_k\right]{\bf P}_k^-\left[{\bf I}-{\bf K}_k{\bf H}_k\right]^T+{\bf K}_k{\bf R}{\bf K}_k^T\)       (8)

여기서, \({\bf K}_k\)는 칼만이득이며, \({\bf H}_k\)는 비선형 측정식 \({\bf h}\left({\bf x}_k\right)\)를 \({\hat{{\bf x}}}_k^-\)로 1차 선형화한 자코비언 행렬이다.

\({\bf H}_k=\left[\left.\frac{\partial{\bf h}\left({\bf x}_k\right)}{\partial{\bf x}_k}\right|_{{\bf x}_k={\hat{{\bf x}}}_k^-}\right]\)       (9)

2.2 Iterated EKF

EKF는 비선형 측정식인 \({\bf h}({\bf x}_k)\)를 테일러 급수로 전개했을 때 1차항만 고려하기 때문에 측정식의 비선형 특성이 매우 크거나 상태변수의 추정오차가 커져서 \({\bf H}_k\)의 선형화 오차가 큰 경우 추정성능이 매우 저하되거나 필터가 발산하는 경우가 발생할 수 있다 (Park &Park 2017). IEKF는 \(k\) 시점에서 동일한 측정치를 이용하여 칼만필터의 측정치 갱신을 반복적으로 수행함으로써 측정모델의 비선형 특성에 의한 자코비언 행렬 계산 오차를 감소시켜 상태변수의 추정성능을 향상시키는 기법이다 (Havlik & Straka 2015).

IEKF의 측정치 갱신을 위해 먼저 상태변수의 초기값을 Eq. (10)과 같이 시간 전파된 상태변수 추정치로 설정한다 (Simon 2006).

\({\hat{{\bf x}}}_k^{(0)}={\hat{{\bf x}}}_k^-\ \)       (10)

측정치 갱신을 반복 수행하여 계산된 \(i\)번째 상태변수 추정치를 \({\hat{\bf x}}_k^{(i)}\)라고 정의하면, \({\hat{\bf x}}_k^{(i)}\)로 선형화된 자코비언 행렬 \({\bf H}_k^{(i)}\)과 칼만이득 \({\bf K}_k^{(i)}\), 상태변수 추정치는 Eqs. (11-13)과 같이 계산된다.

\({\bf H}_k^{(i)}=\left[\left.\frac{\partial{\bf h}\left({\bf x}_k\right)}{\partial{\bf x}_k}\right|_{{\bf x}_k={\hat{{\bf x}}}_k^{(i)}}\right]\)       (11)

\({\bf K}_k^{(i)}={\bf P}_k^-{\bf H}_k^{(i)T}\left[{\bf H}_k^{(i)}{\bf P}_k^-{\bf H}_k^{(i)T}+{\bf R}\right]^{-1}\)       (12)

\({\hat{{\bf x}}}_k^{(i+1)}={\hat{{\bf x}}}_k^-+{\bf K}_k^{(i)}\left[{\bf z}_k-{\bf h}\left({\hat{{\bf x}}}_k^{\left(i\right)}\right)-{\bf H}_k^{(i)}\left({\hat{{\bf x}}}_k^--{\hat{{\bf x}}}_k^{\left(i\right)}\right)\right]\)       (13)

 IEKF의 측정치 갱신은 다음과 같이 연속된 두 개 추정치의 차이가 작거나, 사전 정의된 최대 반복 횟수에 도달할 때까지 Eq. (11)에서 Eq. (13)을 반복적으로 수행한다.

\(xk(i+1)-xk(i)≤ \epsilon\)       (14)

여기서, \(\epsilon\)은 사용자가 정의하는 임계값이다. 마지막 \(i\)번째 반복 추정 결과로부터 최종 측정치 갱신된 IEKF의 상태변수 추정치 및 오차 공분산은 Eqs. (15-16)과 같다.

\({\hat{{\bf x}}}_k^+={\hat{{\bf x}}}_k^{(i)}\)       (15)

\({\bf P}_k^+=\left[{\bf I}-{\bf K}_k^{(i)}{\bf H}_k^{(i)}\right]{\bf P}_k^-\left[{\bf I}-{\bf K}_k^{(i)}{\bf H}_k^{(i)}\right]^T+{\bf K}_k^{(i)}{\bf R}{\bf K}_k^{(i)T}\)       (16)

2.3 Iterated EKF with Consider State

2.3.1 Consider Kalman filter

CKF는 가관측성이 낮은 상태변수를 추정하지는 않지만 오차특성을 공분산 계산에 반영하는 기법으로 상태변수 추정치와 칼만필터 공분산이 유사하여 필터 일관성을 확보할 수 있다 (Hough 2011). 본 연구에서는 측정치에 바이어스 형태의 오차가 포함되는 경우를 가정하여 Eq. (1)에서 상태변수 \( {\bf x}\)를 Eq. (17)과 같이 정의한다.

\({\bf x}=\left[\begin{matrix}{\bf x}_f^T&{\bf x}_c^T\\\end{matrix}\right]^T\)       (17)

여기서, \({\bf x}_f\)는 추정하고자 하는 상태변수이고 \({\bf x}_c\)는 오차 특성만 반영하고자 하는 상태변수이다. 이로부터 \({\bf P}\)와 \({\bf H}\)는 Eq. (18)과 같이 표현할 수 있다.

\({\bf P}=\left[\begin{matrix}{\bf P}_{ff}&{\bf P}_{fc}\\{\bf P}_{cf}&{\bf P}_{cc}\\\end{matrix}\right]\ ,\ \ \ \ \ \ \ \ {\bf H}=\left[\begin{matrix}{\bf H}_f&{\bf H}_c\\\end{matrix}\right]\)       (18)

측정치 잔차 공분산을 \({\bf C}={\bf H}{\bf P}^-{\bf H}^T+{\bf R}\)로 정의하고 Eq. (6)의 칼만이득과 Eq. (8)의 측정치 갱신된 오차 공분산을 분할된 형태로 표현하면 Eq. (19)와 같다 (Zanetti & D’Souza 2013). 편의상 아래첨자 k는 생략하였다.

\({\bf K}=\left[\begin{matrix}{\bf K}_f\\{\bf K}_c\\\end{matrix}\right]=\left[\begin{matrix}\left({\bf P}_{ff}^-{\bf H}_f^T+{\bf P}_{fc}^-{\bf H}_c^T\right){\bf C}^{-1}\\\left({\bf P}_{cf}^-{\bf H}_f^T+{\bf P}_{cc}^-{\bf H}_c^T\right){\bf C}^{-1}\\\end{matrix}\right]\)       (19)

\({\bf P}_{ff}^+={\bf P}_{ff}^--{\bf K}_f{\bf H}\left[\begin{matrix}{\bf P}_{ff}^-\\{\bf P}_{cf}^-\\\end{matrix}\right]-\left[\begin{matrix}{\bf P}_{ff}^-\\{\bf P}_{cf}^-\\\end{matrix}\right]^T{\bf H}^T{\bf K}_f^T+{\bf K}_f{\bf C}{\bf K}_f^T\)       (20)

\({\bf P}_{fc}^+={\bf P}_{fc}^--{\bf K}_f{\bf H}\left[\begin{matrix}{\bf P}_{fc}^-\\{\bf P}_{cc}^-\\\end{matrix}\right]-\left[\begin{matrix}{\bf P}_{ff}^-\\{\bf P}_{cf}^-\\\end{matrix}\right]^T{\bf H}^T{\bf K}_c^T+{\bf K}_f{\bf C}{\bf K}_c^T\)       (21)

\({\bf P}_{cc}^+={\bf P}_{cc}^--{\bf K}_c{\bf H}\left[\begin{matrix}{\bf P}_{fc}^-\\{\bf P}_{cc}^-\\\end{matrix}\right]-\left[\begin{matrix}{\bf P}_{fc}^-\\{\bf P}_{cc}^-\\\end{matrix}\right]^T{\bf H}^T{\bf K}_c^T+{\bf K}_c{\bf C}{\bf K}_c^T\)       (22)

Eq. (19)의 \({\bf K}_f\)를 Eqs. (20-21)에 대입하면,

\({\bf P}_{ff}^+={\bf P}_{ff}^--{\bf K}_f{\bf C}{\bf K}_f^T\)       (23)

\({\bf P}_{fc}^+={\bf P}_{fc}^--{\bf K}_f{\bf H}\left[\begin{matrix}{\bf P}_{fc}^-\\{\bf P}_{cc}^-\\\end{matrix}\right]\)       (24)

Eqs. (23-24)를 살펴보면 \(x_f\)의 공분산 \(P_{ff}^+\)와 \(x_f\)와 사이의 공분산 \(P_{fc}^+\)가 \(K_c\)값에 관계없이 결정됨을 볼 수 있다. 즉, \(K_c\)를 0으로 설정하더라도 \(P_{ff}^+\)와 \(P_{fc}^+\)에는 \(x_c\)의 오차특성이 반영된다. 결과적으로 \(x_f\)의 실제 오차 추정특성이 \(P_{ff}\)에 적절히 반영되어 필터의 일관성이 유지됨을 알 수 있다.

CKF는 Eq. (19)에서 \(K_c=0\)으로 설정하여 \(P_{cc}^+ = P_{cc}^-\)로 유지하고 나머지 오차 공분산은 Eqs. (23)과 (24)를 통해 측정치 갱신을 수행한다. 상태변수에 대한 측정치 갱신은 로 유지하고 에 대한 갱신을 Eq. (25)와 같이 수행하며, 일반적으로 \(x_c\)의 초기값을 알 수 없기 때문에 \({\hat{{\bf x}}}_c^0=0\)이므로 \({\bf h}({\hat{{\bf x}}}_c^-)\)는 생략된다.

\({\hat{{\bf x}}}_f^+={\hat{{\bf x}}}_f^-+{\bf K}_f\left({\bf z}-{\bf h}\left({\hat{{\bf x}}}_f^-\right)-{\bf h}\left({\hat{{\bf x}}}_c^-\right)\right)\)       (25)

2.3.2 Iterated consider EKF

IEKF에 대해서도 가관측성이 낮은 상태변수를 추정하지는 않고 오차특성을 공분산 계산에 반영하여 필터의 일관성을 유지하는 ICEKF를 구현할 수 있다. ICEKF에서 상태변수에 대한 측정치 갱신은 \(x_f\)에 대해서만 수행되므로 Eq. (26)과 같이 초기화 한다.

\({\hat{{\bf x}}}_f^{(0)}={\hat{{\bf x}}}_f^-\)        (26)

반복 수행되는 측정치 갱신식은 Eqs. (27-29)와 같다.

\({\bf H}^{(i)}=\left[\begin{matrix}{\bf H}_f^{(i)}&{\bf H}_c^{(i)}\\\end{matrix}\right],\ \ \ {\bf H}_f^{(i)}=\left[\left.\frac{\partial{\bf h}\left({\bf x}_f\right)}{\partial{\bf x}_f}\right|_{{\bf x}_f={\hat{{\bf x}}}_f^{(i)}}\right],\ \ \ {\bf H}_c^{(i)}=\left[\left.\frac{\partial{\bf h}\left({\bf x}_c\right)}{\partial{\bf x}_c}\right|_{{\bf x}_c=0}\right]\)       (27)

\({\bf K}^{\left(i\right)}={\bf P}^-{\bf H}^{\left(i\right)T}\left[{\bf H}^{\left(i\right)}{\bf P}^-{\bf H}^{\left(i\right)T}+{\bf R}\right]^{-1}=\left[\begin{matrix}{\bf K}_f^{(i)}\\{\bf K}_c^{(i)}\\\end{matrix}\right]\equiv\left[\begin{matrix}{\bf K}_f^{(i)}\\\mathbf{0}\\\end{matrix}\right]\)       (28)

\({\hat{{\bf x}}}_f^{(i+1)}={\hat{{\bf x}}}_f^-+{\bf K}_f^{(i)}\left[{\bf z}-{\bf h}\left({\hat{{\bf x}}}_f^{\left(i\right)}\right)-{\bf H}_f^{(i)}\left({\hat{{\bf x}}}_f^--{\hat{{\bf x}}}_f^{\left(i\right)}\right)\right]\)       (29)

마지막 \(i\)번째 반복 추정 결과로부터 최종 측정치 갱신된 ICEKF의 상태변수 추정치 및 오차공분산은 Eqs. (30-31)과 같다.

\({\hat{{\bf x}}}_f^+={\hat{{\bf x}}}_f^{(i)}\)       (30)

\({\bf P}^+=\left[{\bf I}-{\bf K}^{\left(i\right)}{\bf H}^{\left(i\right)}\right]{\bf P}^-\left[{\bf I}-{\bf K}^{\left(i\right)}{\bf H}^{\left(i\right)}\right]^T+{\bf K}^{\left(i\right)}{\bf R}{\bf K}^{\left(i\right)T}\)       (31)

3. BEARING-ONLY LOCALIZATION OF GNSS INTERFERENCE

 2장에서 제시된 ICEKF를 방위각 도래각을 이용한 GNSS 전파교란원의 2차원 위치추정에 적용하여 성능을 비교하였다. 교란원은 정지된 상태에서 GNSS 교란신호를 송출하는 것으로 가정하였고 센서는 배열안테나 기반 신호처리를 통해 교란신호의 방위각을 측정하고 여러 센서들이 고정된 위치에 일정 간격으로 이격 배치되어 측정된 다수의 방위각 정보를 이용하여 교란원의 위치를 추정하는 것으로 가정하였다.

3.1 Measurement Model for Bearing-only Localization

교란신호의 방위각 도래각을 Fig. 1과 같이 정북 기준으로 교란원 \({\bf p}=\left[\begin{matrix}p_x&p_y\\\end{matrix}\right]^T\)으로부터 센서 \({\bf s}_i=\left[\begin{matrix}s_{xi}&s_{yi}\\\end{matrix}\right]^T\)에 도달하는 각으로 정의하면 \(i\)번째 센서의 도래각 \(\theta_i\)는 Eq. (32)와 같이 표현할 수 있다.

 

Fig. 1. Definition of bearing measurement.

\(\theta_i=\tan^{-1}{\frac{p_x-s_{xi}}{p_y-s_{yi}}}\equiv h_i\)       (32)

센서에서 측정된 도래각에 바이어스 형태의 오차 \(b_i\)와 잡음오차 \(\eta\)이 포함되는 경우 도래각 측정치 \(z_i\)는 Eq. (33)과 같이 표현된다.

\(z_i=\theta_i+b_i+\eta\)        (33)

도래각 측정치의 바이어스 오차 \(b\)는 Eq. (34)와 같이 상관시간이 \(\tau\)이고 평균제곱 오차가 \(\sigma_b^2\)인 1차 마코프 오차로 모델링 하였다 (Brown & Hwang 1997).

\(\dot{b}=-\frac{1}{\tau}b+w_b\)       (34)

여기서, \(w_b\)의 분산 \(q_b={2\sigma_b^2}{\tau}\)이다.

3.2 Reduced-order EKF and IEKF Model

교란원의 위치추정을 위한 칼만필터의 상태변수로 교란원의 위치와 도래각 측정치에 포함된 바이어스 오차를 포함할 수 있지만, 정지상태 교란원에 대한 측정치의 바이어스 오차는 거의 추정되지 않기 때문에 측정치의 바이어스 오차를 제외하고 교란원의 위치만 상태변수로 정의하는 차수감소 필터를 설계할 수 있다.

상태변수를 \({\bf x}=\left[\begin{matrix}p_x&p_y\\\end{matrix}\right]^T\)로 정의하면, 정지상태 교란원에 대한 위치추정이므로 Eq. (1)의 시스템 모델은 Eq. (35)와 같다.

\(\mathrm{\Phi}=I_{2\times2},\ \ \ \ w=0\)       (35)

그리고 \(i\)번째 센서에 대한 Eq. (2)의 비선형 측정모델과 Eq. (9)의 자코비언 행렬은 Eqs. (36-37)과 같이 표현할 수 있다.

\(h_i\left({\bf x}\right)=\theta_i,\ \ \ \ v_i=b_i+\eta\)       (36)

\({\bf H}_{\bf i}=\left[\begin{matrix}\frac{\partial h_i}{\partial p_x}&\frac{\partial h_i}{\partial p_y}\\\end{matrix}\right]=py-syip-s2-px-sxip-s2\)       (37)

3.3 ICEKF Model with Consider State

측정치의 바이어스 오차에 의한 영향을 칼만필터의 공분산 계산에 반영하는ICEKF의 상태변수는 교란원의 위치와 \(N\)개 센서 측정치의 바이어스 오차를 포함하여 Eq. (38)과 같이 정의된다.

\({\bf x}=\left[\begin{matrix}p_x&p_y&b_1&\cdots&b_N\\\end{matrix}\right]^T\)       (38)

여기서, 추정하고자 하는 상태변수 \({\bf x}_{\bf f}=\left[\begin{matrix}p_x&p_y\\\end{matrix}\right]^T\)는 교란원의 위치이고 오차 특성만 반영하고자 하는 상태변수 \({\bf x}_{\bf c}=\left[\begin{matrix}b_1&\cdots&b_N\\\end{matrix}\right]^T\)는 방위각 측정치의 바이어스 오차이다. 시스템 모델은 Eqs. (34-35)로부터 Eq. (39)와 같이 표현할 수 있다.

\(\mathrm{\Phi}=\left[\begin{matrix}1&0&0&0&0\\0&1&0&0&0\\0&0&1-{1}{\tau}&0&0\\\vdots&\vdots&\vdots&\ddots&\vdots\\0&0&0&0&1-{1}{\tau}\\ \end{matrix}\right]\ ,\ \ \ \ w=\left[\begin{matrix}0\\0\\w_{b1}\\⋮\\w_{bN}\end{matrix}\right]\)       (39)

그리고 \(i\)번째 센서에 대한 Eq. (18)의 \(x_f\)와 \(x_c\)에 대한 자코비언 행렬은 Eq. (40)과 같다.

\(H_{fi}={p_y-s_{yi} \over \Vert p-s\Vert^2}- {p_x-s_{xi} \over \Vert p-s\Vert^2} , H_{ci}=⋯1⋯0]\)       (40)

ICEKF의 측정잡음은 센서의 잡음오차로만 정의되므로 Eq. (41)과 같다.

\(v_i=\eta\)        (41)

3.4 Kalman Filter Initialization

칼만필터로 교란원의 위치를 반복 추정하기 위해서는 상태변수와 오차 공분산의 초기값 설정이 필요하다. 상태변수의 초기값을 임의의 값으로 설정하는 경우 실제 교란원의 위치와 임의값의 차이가 크면 측정식의 비선형특성에 의해 필터가 수렴하지 않을 수도 있기 때문에 상태변수의 초기값을 센서에서 측정한 첫번째 방위각 도래각을 이용하여 설정한다. Fig. 1에서 센서 \(s\)로부터 교란원 \(p\)로의 직선을 Line Of Bearing (LOB)라 하면, 센서 \(s_i\)에서 측정되는 LOB의 단위벡터는 Eq. (42)와 같다 (Grabbe et al. 2013).

\({\bf u}_i=\left[\begin{matrix}\sin{z_i}&\cos{z_i}\\\end{matrix}\right]^T\)       (42)

센서  \(s_i\)와 \(s_j\)로 측정되는 두 개의 LOB는 교란원 \(p\)에서 교점을 형성하기 때문에 Eq. (43)의 관계를 만족한다.

\({\bf p}={\bf s}_i+r_i{\bf u}_i={\bf s}_j+r_j{\bf u}_j\)       (43)

여기서, \(r_i\)와 \(r_j\)는 센서와 교란원 사이의 거리이며, Eq. (43)으로부터 Eq. (44)와 같이 구해진다.

\(\left[\begin{matrix}r_i&r_j\\\end{matrix}\right]^T\ =\left[\begin{matrix}{\bf u}_i&-{\bf u}_j\\\end{matrix}\right]^{-1}\left[\begin{matrix}{\bf s}_j&{\bf s}_i\\\end{matrix}\right]\)          (44)

Eqs. (43-44)로부터 교란원의 위치는 \({\bf p}\left(0\right)={\bf s}_i+r_i{\bf u}_i\)로 계산될 수 있으며, 이 값을 칼만필터 상태변수의 초기값으로 정의한다.

칼만필터의 위치 상태변수에 대한 오차공분산의 초기값은 계산된 상태변수 초기값의 오차 특성을 적절히 반영하여야 한다. 일반적으로 위치탐지 오차는 교란원과 센서들간의 기하학적 배치로 정의되는 DOP와 측정오차의 곱으로 정의되기 때문에 (Dempster 2006) 오차공분산의 초기값을 Eq. (45)와 같이 정의하였다.

\({\bf P}\left(0\right)=\left[{\bf H}^T{\bf H}\right]^{-1}\sigma_\theta^2\)       (45)

여기서, \(H\)는 Eq. (37)의 자코비언 행렬로 앞에서 구해진 상태변수 초기값을 이용하여 계산되며 \(\sigma_\theta\)는 도래각 측정오차이다.

4. SIMULATION AND RESULTS

시뮬레이션 성능 분석을 위해 3개의 센서를 Fig. 2와 같이 \(x\)축상에 2 km 간격으로 배치하고 교란원은 센서의 전방에 존재하는 것으로 가정하였다. 각각의 센서는 GPS L1 주파수의 반파장 간격으로 5개의 안테나 소자가 선형으로 배치된 배열안테나로 방위각 도래각을 측정하는 것으로 가정하였다. 먼저 도래각 기반 위치추정의 오차 특성을 확인하기 위해 교란원의 위치에 따른 타원오차를 계산하였다. Eq. (45)로부터 타원오차의 크기는 공분산행렬 고유값의 제곱근, 방향은 고유벡터로 정의되며, 도래각 측정오차를 0.7도로 가정하고 6개 지점에 대해 각각의 교란원 위치에서의 타원오차를 Fig. 2에 나타 내었다. 각각에 대한 위치 추정오차의 distance root mean square (drms)는 J2가 40.3 m, J3가 154.6 m, J5가 24.8 m, J6가 205.5 m로 센서에서 거리가 멀수록 위치 추정오차가 커지며 J6와 같이 센서들의 기저선과 교란원의 LOB가 이루는 각이 작아질수록 오차가 커지는 것을 알 수 있다. 특히 J6의 경우 위치 추정값이 길쭉한 타원상에 형성되기 때문에 실제 위치와 추정위치의 차이가 크게 발생할 수 있으며 이로 인해 측정식의 비선형특성에 의한 영향을 많이 받게 됨을 유추할 수 있다.

Fig. 2. Error ellipse of bearing-only localization with 3 sensors for each interference location.

차수감소 EKF와 IEKF 및 ICEKF의 성능을 비교하기 위해 몬테칼로 시뮬레이션을 수행하였다. 방위각 도래각 측정치의 1차 마코프 바이어스 오차 \(\sigma_b\)는 0.5도, \(\tau\)는 1000초로 가정하고 잡음오차는 0.5도로 가정하였다. 교란원은 DOP이 좋은 위치인 J2와 DOP이 나쁜 위치인 J6에 있는 것으로 가정하고 각각의 교란원에 대해 몬테칼로 시뮬레이션을 100회 수행하여 위치 추정오차의 drms와 칼만필터에서 계산한 오차 공분산을 비교하였다.

Figs. 3과 4는 각각 교란원이 J2와 J6에 있는 경우에 대해 EKF와 IEKF, ICEKF의 위치 추정 결과를 나타낸 것이고 1은 마지막 시점에서의 3개 필터의 추정오차와 칼만필터 공분산을 비교한 것이다. 위치 추정오차 관점에서 교란원이 J2에 있는 경우에는 추정오차가 크지 않아 측정식의 비선형특성에 의한 영향이 작기 때문에 3개 필터의 성능이 거의 유사한 반면, 교란원이 J6에 있는 경우 큰 추정오차로 인해 측정식의 비선형특성 영향이 커져 EKF는 발산하는 경우가 생기지만 IEKF와 ICEKF는 발산하지 않고 교란원의 위치를 추정함을 볼 수 있다. 한편, 위치 추정오차와 칼만필터 공분산에 대한 필터의 일관성 관점에서 살펴보면 측정치의 바이어스 오차를 측정잡음에 반영한 차수감소 EKF와 IEKF는 실제 추정오차에 비해 칼만필터 공분산이 더 작게 계산되어 필터의 일관성을 확보하지 못하는 반면, ICEKF는 칼만필터에서 계산된 공분산이 실제 추정오차와 거의 유사한 값으로 계산되어 필터의 일관성이 확보되는 것을 확인할 수 있다.

Fig. 3. Localization error for J2 (Gray line: each 100 runs, thick line: drms error, dotted line: STD of Kalman filter).

Fig. 4. Localization error for interference at J6 (Gray line: each 100 runs, thick line: drms error, dotted line: STD of Kalman filter).

1. Comparison of drms error and STD of Kalman filter for interference J2 and J6.

Filter	Interference at J2		Interference at J6
	drms	STD of Kalman filter	drms	STD of Kalman filter
EKF (m)  IEKF (m)  ICEKF (m)	30.8 30.8 30.5	8.0 8.0 27.3	12021.1 152.2 142.7	44.1 41.1 135.4

 

IEKF와 ICEKF의 실제 추정오차를 비교해 보면 거의 유사한 것을 볼 수 있다. 이는 칼만필터에서 연속적인 측정치 갱신을 통해 잡음오차에 의한 영향은 필터링되어 감소되지만, 바이어스 오차에 의한 영향은 그대로 남아 위치 추정오차에 반영되기 때문이다. 반면, IEKF와 ICEKF의 공분산을 비교해 보면 IEKF의 공분산은 실제 추정오차 보다 훨씬 작게 계산된다. 이는 IEKF에서 측정치의 바이어스 오차가 측정잡음에 포함됨으로써 필터 관점에서는 측정잡음이 필터링되어 위치 추정오차가 작아진 것으로 계산되었기 때문이다. ICEKF는 바이어스 오차의 공분산이 위치 오차 공분산 계산에 반영되기 때문에 실제 추정오차와 필터의 공분산이 유사하게 계산된다.

5. CONCLUSIONS

본 논문에서는 방위각 도래각을 이용하여 GNSS 교란원의 위치를 추정하는 기법으로 도래각 측정식의 비선형특성과 측정치에 포함된 바이어스 오차의 특성을 반영하는 ICEKF를 설계하고 시뮬레이션을 통해 EKF, IEKF와 성능을 비교하였다. CKF는 가관측성이 낮거나 시스템 성능에 영향을 미치지 않는 상태변수를 추정하지 않고 그 오차특성만 필터의 공분산 계산에 반영하는 필터이다. 그리고 IEKF는 상태변수 초기값의 큰 오차에 의해 비선형 도래각 측정치의 선형화 오차가 커져 EKF가 발산하는 것을 방지하기 위해 동일한 측정치로 측정치 갱신을 반복하는 기법이다. ICEKF는 CKF를 IEKF에 적용한 것으로 선형화 오차에 의한 필터의 발산 억제와 추정 성능 향상 및 추정치와 필터 공분산이 일치하는 필터의 일관성을 확보하는 장점을 가진다. 지상기반 2차원 위치추정에 대한 시뮬레이션을 수행하여 도래각 측정치에 포함된 바이어스 오차의 특성을 필터 공분산 계산에 반영한 ICEKF의 성능이 EKF 및 IEKF에 비해 향상됨을 확인하였다.

AUTHOR CONTRIBUTIONS

Conceptualization, Y. Park and K. Song; methodology, Y. Park and K. Song; validation, Y. Park; writing-original draft preparation, Y. Park; writing-review and editing, K. Song.

CONFLICTS OF INTEREST

The authors declare no conflict of interest.

BIOGRAPHY

Youngbum Park received the Doctor’s degree in Aerospace engineering from Seoul National University in 2017. Since 2001, he has been working for Agency for defense development. His research interests include estimation, localization, and integrated navigation system.

Kiwon Song received the Doctor’s degree in Electronics engineering from Chung-nam National University in 2002. His research interests include design of satellite navigation system architecture, GNSS signal generation processing, and integration optimal filter of INS/GPS.