Along-Track Position Error Bound Estimation using Kalman Filter-Based RAIM for UAV Geofencing

Abstract

Geofencing supports unmanned aerial vehicle (UAV) operation by defining stay-in and stay-out regions. National Aeronautics and Space Administration (NASA) has developed a prototype of the geofencing function, SAFEGUARD, which prevents stayout region violation by utilizing position estimates. Thus, SAFEGUARD depends on navigation system performance, and the safety risk associated with the navigation system uncertainty should be considered. This study presents a methodology to compute the safety risk assessment-based along-track position error bound under nominal and Global Navigation Satellite Systems (GNSS) failure conditions. A Kalman filter system using pseudorange measurements as well as pseudorange rate measurements is considered for determining the position uncertainty induced by velocity uncertainty. The worst case pseudorange and pseudorange rate fault-based position error bound under the GNSS failure condition are derived by applying a Receiver Autonomous Integrity Monitor (RAIM). Position error bound simulations are also conducted for different GNSS fault hypotheses and constellation conditions with a GNSS/INS integrated navigation system. The results show that the proposed along-track position error bounds depend on satellite geometries caused by UAV attitude change and are reduced to about 40% of those of the single constellation case when using the dual constellation.

1. INTRODUCTION

Geo-fencing은 무인항공기의 안전성을 위협할 수 있는 문제를 최소화하기 위해 전체 공역을 기체가 운용을 수행하는 영역인 Stay-in 영역과 침범해서는 안 되는 영역인 Stay-out 영역으로 분할하여 운용하는 방식을 말한다. NASA에서는 geo-fencing 운용에서 Stay-out 영역으로의 침범을 방지하여 운용의 안전성을 보장하기 위해 geo-fencing function의 프로토타입인 SAFEGUARD를 개발하였다 (Dill et al. 2016). SAFEGUARD는 무인항공기의 위치 추정치와 geo-fencing 영역의 경계를 비교하여 Stay-out 영역으로의 침범을 방지하기 때문에 항법 시스템 성능에 크게 의존한다. 사전 연구에서는 항법 불확실성으로 인해 발생할 수 있는 SAFEGUARD 기반 geo-fencing 운용의 위험상황을 방지하기 위해 위험성 확률 평가 기반 위치 오차 경계를 사용하여 Stay-out 영역으로의 침범을 방지하는 운용 방식이 제안되었다 (Kim et al. 2019). 하지만 Kim et al. (2019)에서 제안된 방식을 기반으로 무인항공기의 geo-fencing 운용 안전성을 보장하는 데에는 다음 두 가지 한계점이 있다: 1) 항법 시스템 정상상태만을 고려한 위험성 확률 평가를 수행하였으며, 2) GPS 단일 항법 시스템을 가정하여 환경변화에 강건한 항법을 위해 복합 항법 시스템을 사용하는 일반적인 무인항공기 항법 시스템에 대한 고려가 부족하다.

따라서 본 연구에서는 사전연구에서 제안된 geo-fencing 운용의 위치 오차 경계 산출 과정에서 항법 시스템의 고장상태를 고려할 수 있는 방법론을 제안하였으며, 복합 항법 시스템 사용 시 발생할 수 있는 속도 불확실성으로 인한 위치 오차를 고려한 기체의 속도 방향 (종뱡향, Along-Track) 위치 오차 경계 산출 방법론을 제안하였다. 속도 불확실성에 의한 종방향 위치 오차를 고려하기 위해서는 위치와 속도 오차 간의 연관성을 고려해야한다. 따라서 본 연구에서는 속도와 위치를 동시 추정하고, 속도 추정치를 기반으로 위치를 업데이트하는 칼만필터 기반의 종방향 위치 오차 모델을 활용하였다. 또한, 항법 시스템 고장에 대한 안전성 보장을 위해 칼만필터 이노베이션을 검정 통계량으로 사용하는 Receiver Autonomous Integrity Monitor (RAIM)를 적용하였으며, 정확한 속도 산출을 위해 의사거리 측정치와 함께 의사거리 변화율 측정치를 사용하는 경우에 대한 위치 오차 경계 산출 방법론을 제안하였다.

2장에서는 위험성 확률 평가 기반 위치 오차 경계 산출 방법, 3장에서는 속도 불확실성을 고려한 위치 오차 경계 산출 방법 및 시뮬레이션 결과 그리고 4장에서는 본 논문의 요약 및 결론을 기술한다.

본 논문은 2022 항법시스템학회 정기학술대회에서 발표한 Nam et al. (2022)의 내용을 확장하여 작성하였다.

2. SAFETY RISK ASSESSMENT BASED POSITION ERROR BOUND

2.1 Concept of SAFEGUARD

Geo-fencing은 무인 항공기 운용의 안전성을 위협할 수 있는 문제를 최소화하기 위해, 전체 공역을 크게 둘로 분할하여 운용하는 방식을 말한다. 하나는 기체가 운용을 수행하는 영역으로 Stay-in 영역, 다른 하나는 기체가 침범해서는 안 되는 영역으로 Stay-out 영역이다. 미항공우주국 (NASA)에서는 geo-fencing 운용의 안전성을 보장하기 위해 Stay-out 영역으로의 침범을 방지할 수 있는 geo-fencing function의 프로토타입인 SAFEGUARD를 개발하였다. SAFEGUARD는 Fig. 1과 같이 Stay-out 영역으로의 침범을 효과적으로 방지하기 위해 Stay-in 영역과 Stay-out 영역 사이에 세 가지 한계선을 정의하여 사용한다 (Dill et al. 2016).

Fig. 1. Boundaries and regions in SAFEGUARD (left) and boundaries for geospatial constraint conformance (right) (Dill et al. 2016).

세 한계선은 Stay-out 영역에서부터 절대 침범해서는 안되는 최종 한계선 (Hard boundary), 최종 한계선의 침범을 막기 위해 무인항공기를 강제 추락시키는 추락 한계선 (Terminate boundary) 그리고 Stay-out 영역에 침범할 가능성이 있는 경우 사용자에게 경고를 보내는 경고 한계선 (Warning boundary)이다. 최종 한계선은 Stay-out 영역으로부터 결정되며 최종 한계선을 기준으로 기체의 운동상태를 고려하여 추락 한계선의 크기(\(\varepsilon\))가 운용 중에 정해진다. 마지막으로 경고 한계선의 크기(\(\rho\varepsilon\))는 추락 한계선을 기준으로 기체의 경고 신호에 대한 반응 속도를 고려하여 결정된다. SAFEGUARD는 각 한계선과 사용자의 위치 추정치를 비교하여 접촉한 경우 각 한계선의 역할을 수행하여 Stay-out 영역으로의 침범을 방지한다.

2.2 Safety Risk Assessment

SAFEGUARD 기반 운용의 안전성은 항법 시스템 성능에 크게 의존하기 때문에, 항법 시스템의 불확실성을 반드시 고려해야 한다. Kim et al. (2019)에서는 SAFEGUARD와 같은 geo-fencing function을 사용한 무인항공기의 항법 불확실성을 고려하기 위해 항법 시스템으로부터 산출된 위치해가 아닌 위치 오차 경계 기반의 geo-fencing 운용 방식을 제안하였다. 위치 오차 경계를 산출하기 위해 운용 중 발생할 수 있는 인명피해 확률을 고려하는 위험성 확률 평가 기반 위치 오차 경계 산출 방법론을 제안하였다. 구체적인 임무 경로를 고려하지 않고 경고 한계선 근처에서 운용한다는 가정 하에 비행 기술 오차를 제외한 항법 시스템 오차만을 고려하여 geo-fencing 운용의 위험성 확률을 평가하였으며, 해당 확률이 항법 시스템의 안전성 요구조건을 만족하는 위치 오차 경계를 산출하였다.

Kim et al. (2019)에서 정의한 geo-fencing 운용 중 항법 불확실성으로 인명피해가 발생할 수 있는 상황은 Fig. 2와 같다. Fig. 2의 붉은 원은 기체의 실제 위치, 파란 원은 기체의 추정 위치, 그리고 파란색 박스는 위치 오차 경계를 나타낸다. Case 1은 실제 위치가 경고 한계선과 최종 한계선 사이에 위치해 있으면서 경고 신호에 반응하지 않고 추락 한계선에 의해 Stay-out 영역 내부에 추락하는 경우, Case 2는 실제 위치가 최종 한계선 내에 위치한 경우, 그리고 Case 3는 위치 오차 경계 내에 실제 위치가 존재하지만 경고 한계선의 신호에 반응하지 않아 추락 한계선에 의해 지상에 충돌하는 경우를 나타낸다. Geo-fencing 운용의 위험성 확률은 세 경우의 확률의 합으로 계산되며, 세 가지 추가적인 가정을 통해 해당 확률은 다음과 같이 표현된다 (Kim et al. 2019): 1) 위치 오차 경계가 높은 수준의 안전성 요구조건을 충족한다는 가정하에 위치 오차가 위치 오차 경계를 넘지 않을 확률을 1로 가정, 2) 실제 위치가 최종 한계선을 넘었을 때 인명피해가 발생할 확률을 1로 가정, 그리고 3) 위치 오차가 위치 오차 경계, 한계선들의 간격을 합한 것 보다 클 확률을 0으로 가정하였다.

\(Pr(HH)=Pr(HH|CO_{ground})\cdot Pr⁡(No~React)+Pr⁡(PE>PE~Bound)\cdot Pr⁡(No~React)\)       (1)

여기서 \(Pr(HH)\)는 항법 불확실성에 의한 인명피해 (Harm to Human, HH) 발생 확률, \(Pr(HH|CO_{ground})\)는 무인항공기 지상 충돌\((CO_{ground})\) 시 인명피해 발생 확률, \(Pr⁡(No~React)\)는 무인항공기가 경고 한계선의 신호에 반응하지 않을 확률 그리고 \(Pr⁡(PE>PE~Bound)\)는 위치 오차 (Position Error, PE)가 위치 오차 경계 (PE Bound)보다 클 확률을 나타낸다. Kim et al. (2019)에서는 식 (1)을 기반으로 GPS 단일 항법 시스템의 정상상태 오차모델을 사용하여 인명피해 발생 확률이 항법 시스템의 안전성 요구조건을 만족하는 위치 오차 경계를 산출하였다. 하지만 항법 시스템의 불확실성으로부터 운용의 안전성을 보장하기 위해서는 항법 시스템 고장을 고려한 위험성 확률 평가도 함께 수행되어야 한다. 따라서 본 연구에서는 식 (1)을 기반으로 항법 시스템의 고장상태에 대한 인명피해 확률 평가 기반 위치 오차 경계 산출 방법론을 제안하였다.

Fig. 2. 3 Fatal cases in geofencing operation (Kim et al. 2019).

항법 시스템의 정상상태와 고장상태는 서로 상호 배타적인 사건이다. 따라서 인명피해 발생확률은 항법 시스템의 상태에 따라서 식 (2)와 같이 표현될 수 있다.

\(Pr(HH)=Pr⁡(HH|FF)  Pr⁡(FF)+Pr⁡(HH|F)  Pr⁡(F)\)                                               (2)

여기서 \(Pr⁡(HH|FF)\)는 항법 시스템 정상상태 (Fault Free, FF) 하에서의 인명피해 발생확률, \(Pr⁡(HH|F)\)는 고장상태 (Fault, F) 하에서의 인명피해 발생확률, Pr⁡(FF)는 항법 시스템 정상상태의 사전 발생확률 그리고 \(Pr⁡(F)\)는 항법 시스템 고장의 사전 발생확률을 나타낸다.

식 (1) 우항의 위치 오차가 위치 오차 경계를 넘을 확률은 항법 시스템 상태 따라 다르게 정의되기 때문에 식 (3)과 같이 표현할 수 있다.

\(Pr⁡(HH|FF)  Pr⁡(FF)+Pr⁡(HH|F)  Pr⁡(F) \\ =Pr(HH|CO_{ground} )\cdot Pr⁡(No~React) +Pr⁡(PE>PE Bound|FF)\cdot Pr⁡(FF) \cdot Pr⁡(No~React) \\ +Pr⁡(PE>PE~Bound|F)\cdot Pr⁡(F) \cdot Pr⁡(No~React)\)    (3)

항법 시스템 정상상태 및 고장상태의 사전 발생확률의 합은 항상 1이기 때문에 식 (3)은 식 (4)와 같이 정리할 수 있다.

\(Pr⁡(HH|FF)  Pr⁡(FF)+Pr⁡(HH|F)  Pr⁡(F) \\ =(Pr(HH|CO_{ground} )+Pr⁡(PE>PE~Bound|FF) )\cdot Pr⁡(FF)\cdot Pr⁡(No~React) \\ +(Pr(HH|CO_{ground} )+Pr⁡(PE>PE~Bound|F) )\cdot Pr⁡(F)\cdot Pr⁡(No~React)\)    (4)

따라서 항법 시스템 정상상태 및 고장상태 하에서의 인명피해 발생 확률은 식 (5, 6)을 통해 평가될 수 있다.

\(Pr⁡(HH|FF)=(Pr(HH|CO_{ground} )+Pr⁡(PE>PE~Bound|FF) ) \cdot Pr⁡(No~React)\)                 (5)

\(Pr⁡(HH|F)=(Pr(HH|CO_{ground} )+Pr⁡(PE>PE~Bound|F) )\cdot Pr⁡(No~React)\)                         (6)

식 (5, 6)을 통해 평가된 인명피해 발생 확률은 항법 시스템의 정상상태와 고장상태에 할당된 시스템 전체의 안전성 요구조건 (Target Level of Safety, TLS)을 만족해야 운용의 안전성을 보장할 수 있다. 따라서 본 연구에서는 할당된 안전성 요구조건으로부터 항법 시스템의 무결성 요구조건(\(Pr⁡(PE>PE~Bound)\))을 도출하였으며, 정상상태와 고장상태의 오차모델을 기반으로 무결성 요구조건을 만족하는 위치 오차 경계(보호수준)를 산출하였다. 다음 장에서는 항법 시스템 정상상태와 고장상태에 대한 위치 오차 경계 산출 방법론을 자세하게 기술한다.

3. KALMAN FILTER BASED ALONG-TRACK POSITION ERROR BOUND

Kim et al. (2019)에서 제안된 위치 오차 경계는 속도 불확실성으로 발생할 수 있는 기체 이동경로의 종방향 (Along-track) 위치 불확실성을 고려하지 않는다. 하지만 일반적인 무인 항공기의 항법 시스템은 속도와 위치를 동시에 추정하는 복합 항법 시스템을 사용하기 때문에 속도 불확실성으로 발생할 수 있는 위치 불확실성을 고려해야한다. 따라서 본 연구에서는 칼만필터 기반 GNSS/INS 융합 항법 시스템의 속도 불확실성을 고려한 종방향 위치 오차 경계 산출 방법을 제안하였으며, 고장상태의 경우 GNSS 고장상황만을 고려하였다. GNSS/INS 전체 시스템의 안전성 보장을 위한 INS 고장상황에 대해서는 후속 연구에서 다룰 예정이다.

3.1 Position Error Bound Under Normal Condition

속도 불확실성에 의한 위치 불확실성을 고려하기 위해서는 속도와 위치 간의 연관성을 고려한 위치 오차 모델을 사용해야한다. 따라서 본 연구에서는 속도와 위치를 동시 추정하고, 추정한 속도를 기반으로 위치를 업데이트하는 칼만필터 기반 GNSS/INS 융합 항법 시스템의 위치 오차 모델을 사용하였다. 칼만필터의 시스템 업데이트는 식 (7, 8)과 같다.

\(\hat{x}_k^-=F_{k-1} \hat{x}_{k-1}^+\)                                                                                               (7)

\(P_k^-=F_{k-1} P_{k-1}^+ F_{k-1}^T+Q_{k-1}\)                                                                                 (8)

여기서 \(\hat{x}_k^-\)는 \(k\)번째 에폭의 시스템 업데이트된 추정치로 위치, 속도, GNSS 수신기 시계 오차와 관성센서 오차로 구성된다. \(F_{k-1}\)는 시스템 업데이트 행렬, \(\hat{x}_{k-1}^+\)는 이전 에폭의 추정치 벡터, \(P_k^-\)는 \(\hat{x}_k^-\)의 공분산 행렬, \(P_{k-1}^+\)는 \(\hat{x}_{k-1}^+\)의 공분산 행렬 그리고 \(Q_{k-1}\)는 시스템 업데이트 잡음의 공분산 행렬을 나타낸다. 시스템 업데이트는 추정된 속도를 기반으로 위치를 업데이트하기 때문에 시스템 업데이트된 공분산 행렬을 통해 속도 불확실성으로 인한 위치 오차 불확실성을 고려할 수 있다.

GNSS 측정치 모델과 칼만필터의 측정치 업데이트는 식 (9-11)과 같다. GNSS 측정치의 경우 정확한 속도 추정을 위해 의사거리 측정치와 의사거리 변화율 측정치 함께 사용하는 경우를 가정하였다.

\(z_k=\left[\begin{array}{c}\rho_k \\ \dot{\rho}_k \end{array}\right]=\left[\begin{array}{c}H_{\rho,k} & 0 \\ 0&H_{\dot{\rho},k} \end{array}\right] x_k+v_k=H_k x_k+v_k\)                                                     (9)

\(\hat{x}_k^+=\hat{x}_k^-+K_k (z_k-Η_k \hat{x}_k^- )\)                                                                                 (10)

\(P_k^+=(I-K_k H_k ) P_k^-\)                                                                                     (11)

식 (9)에서 \(z_k\)는 측정치 벡터이며, 가시 위성의 의사거리 측정치(\(\rho_k\))와 의사거리 변화율 측정치(\({\dot{\rho}_k}\))로 구성된다. \(H_{\rho,k}\)와 \(H_{\dot{\rho},k}\)는 각각 의사거리와 의사거리 변화율에 대한 측정치 행렬, \(x_k\)는 상태 벡터, \(v_k\)는 측정치 잡음을 나타낸다. 식 (10, 11)은 추정치와 공분산에 대한 칼만필터 측정치 업데이트를 나타내며, \(K_k\)는 칼만 게인, \(P_k^+\)는 \(\hat{x}_k^+\)의 공분산 행렬을 나타낸다. 본 연구에서 가정한 GNSS/INS 항법 시스템은 의사거리 변화율을 활용해 속도 추정치를 보정하기 때문에 칼만필터의 시스템 업데이트를 통해 의사거리 변화율 측정치의 오차가 위치 오차에 전파된다.

고장이 없는 정상상태의 경우 종방향 위치 오차를 영평균 정규분포로 가정할 수 있으며, 종방향 위치 오차 경계는 식 (12)와 같이 계산된다.

\(PEB_{Nominal}=K_{ff} \cdot\sigma_{along-track}\)                                                                     (12)

여기서 \(K_{ff}\)는 항법 시스템 정상상태에 할당된 안전성 요구조건을 바탕으로 식 (5)를 통해 계산된 정상상태 항법 시스템의 무결성 요구조건에 대응하는 백분위수, \(\sigma_{along-track}\)는 칼만필터의 종방향 위치 추정치의 표준편차를 나타낸다.

3.2 Position Error Bound Under GNSS Fault Condition

칼만필터 기반 GNSS/INS 항법 시스템의 GNSS 고장에 대한 안전성 보장 방법으로 칼만필터의 이노베이션 벡터를 검정 통계량으로 사용하는 RAIM (Tanıl et al. 2018)과 칼만필터의 추정치를 기반으로 계산한 잔차 (Residual) 벡터를 검정 통계량으로 사용하는 RAIM (Bhattacharyya & Gebre-Egziabher 2015)이 제안되었다. 하지만 Tanıl et al. (2018)의 경우 의사거리 측정치 사용과 스푸핑 상황만을 가정하였으며, Bhattacharyya & Gebre-Egziabher (2015)의 경우 의사거리와 의사거리 변화율 측정치의 사용을 모두 고려하였지만 의사거리 측정치 고장이 시간에 따라 선형 함수를 따른다는 가정을 통해 의사거리 변화율 측정치의 고장을 고려하였다. 따라서 본 연구에서는 의사거리 변화율 고장의 영향에 의한 항법 불확실성의 일반적인 고려를 위해 의사거리와 의사거리 변화율 측정치의 고장을 동시에 고려한 최악 상황 고장 (Worst Case Fault) 하의 위치 오차 경계 산출 방법론을 제안한다.

본 연구에서는 이노베이션 벡터를 검정 통계량으로 사용하는 Tanıl et al. (2018)의 RAIM을 GNSS 고장 모니터로 사용하였다. 해당 모니터의 \(k\)번째 에폭의 검정 통계량(\(q_k\))은 식 (13)과 같이 정의된다.

\(q_k=\sum\limits_{j=1}^k \gamma_j^T S_j^{-1} \gamma_j\)                                                                                     (13)

여기서 \(\gamma_j\)는 \(j\)번째 에폭의 이노베이션 벡터 그리고 \(S_j\)는 이노베이션 벡터의 공분산 행렬을 나타낸다. 본 연구에서는 의사거리와 의사거리 변화율 측정치를 함께 고려한 이노베이션 벡터를 사용하였다. 정상상태일 경우 식 (13)은 카이제곱분포를 따르며, 누적된 에폭을 \(k\), 에폭당 측정치 수가 \(M\)이면 해당 카이제곱분포는 \(k\cdot M-5\)의 자유도를 갖는다. 반면, 고장상태일 경우 카이제곱분포는 이노베이션의 바이어스로 인해 non-centrality를 갖게 된다.

최악 상황 고장은 고장에 대한 추가적인 가정 없이 가장 보수적인 상황을 가정하여 고장의 영향을 고려하는 방법론이다. 최악 상황 고장은 미검출 고장으로 인한 위치 오차가 운용의 위치 오차 요구조건을 넘을 확률인 무결성 위협 확률을 최대화하는 고장을 말하며, 해당 확률은 식 (14)와 같이 정의된다.

\(P(|PE|>AL,q_k< t_k| GNSS~Failure)\)                                                             (14)

여기서 Alert Limit (\(AL\))은 경고 수준, \(T_k\)는 모니터의 임계값을 나타낸다. 식 (14)의 확률값은 발생한 GNSS 고장에 따라 달라지며, 위치 오차에 큰 영향을 주면서 검정 통계량에 영향이 미미한 고장일수록 무결성 위협 확률이 증가하는 것을 확인할 수 있다.

본 연구에서는 Tanıl et al. (2018)과 Joerger & Pervan (2013)을 기반으로 GNSS 고장을 바이어스 벡터로 가정하고, 고장 벡터의 영향을 정량화하는 기울기(\(\alpha_k^2\))를 정의하여 이를 최대화하는 고장 벡터의 방향을 결정하였다.

\(\alpha_k^2={(E[\tau_a \tilde{x}_k^+ ])^2 \over \lambda_k^2}={E[\tau_a \tilde{x}_k^+ ])^2 \over \sum\limits_{i=1}^k E[\gamma_i^T ] S_i^{-1} E[\gamma_i ] }\)                                                                           (15)

식 (15)의 분자는 종방향 추정치 오차의 기대값의 제곱이며, \(\tau_a\)는 종방향 오차 성분을 추출하는 벡터, \(\tilde{x}_k^+\)는 \(k\)번째 에폭의 추정치 오차이다. 분모는 검정 통계량의 non-centrality의 제곱이다. 고장 벡터가 위치 오차에 큰 영향을 주면서 검정 통계량에 미치는 영향이 작을수록 기울기가 커진다. 따라서 기울기를 최대화하는 고장 벡터가 최악 상황 고장 벡터가 된다.

고장 벡터가 포함된 측정치 모델을 활용해 식 (15)를 고장 벡터의 함수로 나타낼 수 있으며, 측정치 모델은 식 (16)과 같다.

\(z_k=H_k x_k+v_k+f_k\)                                                                                     (16)

식 (16)은 정상상태의 측정치 모델에 고장 벡터(\(f_k\))가 추가된 형태이다. 이때 식 (16)의 고장 벡터는 위성의 고장 가설, 고장 방향, 그리고 고장의 크기에 따라 다르게 정의된다. 식 (17)은 단일 위성 고장에 대한 고장 벡터를 나타내며, 위성 고장 발생 시 의사 거리와 의사 거리 변화율에 모두 바이어스가 발생한다고 가정하였다.

\(f_k=\left[\begin{array}{c}b_{k,1}\\ ⋮\\ 0\\ \dot{b}_{k,1}\\ ⋮\\ 0\end{array}\right]=\left[\begin{array}{cc}1&0\\ ⋮&⋮\\ 0_\infty& 0\\ 0&1\\ ⋮&⋮ \\ 0& 0 \end{array}\right]\left[\begin{array}{c}b_{1,k}\\ \dot{b}_{1,k} \end{array}\right]=T_1 f_{1,k}\)                                                             (17)

식 (17)은 1번 위성 고장가설에 대한 측정치 고장 벡터를 나타내며, \(b_{1,k}\)은 의사 거리의 바이어스, \(\dot{b}_{1,k}\)은 의사 거리 변화율의 바이어스를 나타낸다. 가시 위성 중 고장 위성 조합만큼 가설이 존재하며, 가설 중에서 가장 큰 최대 기울기를 갖는 가설이 최악 고장 가설이 된다. 고장 가설 하에서 최대 기울기를 구하는 과정은 APPENDIX A에 자세하게 기술되어 있다.

GNSS 고장 하에서 종방향 위치 오차 분포는 정상상태와 달리 영평균을 가정할 수 없으며, 본 논문에서는 최악 고장 상황을 가정해 종방향 위치 오차의 평균을 보수적으로 설정하였다. GNSS 고장 하의 종방향 위치 오차 경계는 식 (18)과 같이 계산된다 (Walter & Enge 1995).

\(PEB_{GNSS~Failure}=\left( \overset{\max}{i}⁡\alpha_{i,k,max}\right ) T_k+K_{gnssf}\cdot \sigma_{along-track}\)                                           (18)

식 (18)에서 \(\alpha_{i,k,max}\)는 \(i\)번째 고장가설의 기울기의 최댓값, \(T_k\)는 이노베이션 모니터의 임계값, \(K_{gnssf}\)는 GNSS 고장상태의 안전성 요구조건과 식 (6)을 통해 계산된 GNSS 고장상태 무결성 요구조건에 대응하는 백분위수를 의미한다.

3.3 Simulation Results

이번 절에서는 앞서 기술된 방법을 기반으로 정상상태와 GNSS 고장상태 하에서 종방향 위치 오차 경계 산출 시뮬레이션을 수행하였다. 식 (1)의 Pr⁡(\(HH|CO_{ground}\))는 대한민국 소도시 인구밀도인 약 70명/km²를 가정하여 Kim et al. (2019)과 동일하게 고려하였다. 항법 시스템에 할당된 TLS는 Grimsley (2004)에서 제시한 \(10^{-6} fatality/hr\)를 가정하였으며, 그 중 종방향에 \(5\cdot 10^{-7}\)을 할당하였다. INS로는 Consumer Grade IMU의 사용을 가정하였으며 (Groves 2013), GNSS로 단일 위성군(GPS)과 다중 위성군(GPS+Galileo)을 활용하였다. GNSS 고장의 사전 발생확률은 도심 운용을 고려하여 US Department of Defense (2020)에서 제시하는 \(10^{-4}/hr\)보다 보수적인 \(10^{-3}/hr\)를 사용하였다. GNSS 고장 가설의 경우 사전 발생확률을 기반으로 이중 위성 고장 가설까지 고려하였다 (Joerger et al. 2014). 그 외의 시뮬레이션 변수는 Table 1에 나타냈으며, 영국 텔스콤 인근(북위 51°, 동경0°)에서 Fig. 3의 비행경로와 속도 프로파일을 가정하여 시뮬레이션을 수행하였다.

Table 1. Simulation parameters.

Parameters	Values (unit)
\(Pr⁡(F)\) \(Pr⁡(No~React)\) \(TLS_{FF}\) \(TLS_F\) \(P_{Continuity~risk}\) \(\sigma_\rho\) \(\sigma_{\dot\rho}\)	\(10^{-3}/hr\) \(10^{-4}/hr\) \(4.8\cdot 10^{-7}/hr\) \(2\cdot 10^{-8}/hr\) \(10^{-6}/hr\) 1 m 0.1 m/s

 

Fig. 3. Simulation trajectory (left) and velocity profile (right) in ECEF.

Fig. 4는 GPS 단일 위성군 사용 상황에 대한 정상상태 및 GPS 고장상태 하에서의 종방향 위치 오차 산출 시뮬레이션 결과(좌)와 ECEF 위치 오차의 표준편차 추정치(우)를 나타낸다. GPS 단일 위성군에서 정상상태 위치 오차 경계는 평균 2 m 정도에서 비교적 일정하게 유지된다. 반면, GPS 고장상태의 위치 오차 경계는 정상상태보다 항상 크며, 이중 고장의 위치 오차 경계(빨간선)가 단일 고장(파란선)의 경우 보다 항상 큰 것을 확인할 수 있다. GPS 고장상태의 경우 시간이 지날 수록 누적되는 고장의 영향이 증가하여, 기체의 방향이 일정한 구간에서 위치 오차 경계가 증가하는 경향을 보이는 것을 확인할 수 있다. 또한, GPS 고장상태에서는 정상상태에 비해 기체의 방향에 따라 오차 경계가 크게 변화한다. Fig. 4의 좌측 그래프를 보면, 14s에서 46s 사이에 기체가 진행 방향을 바꾸면서 위치 오차 경계가 감소함을 확인할 수 있다 (Fig. 3의 시뮬레이션 경로 참고). 이는 기체의 회전으로 기체 기준 좌표계에 대한 위상기하가 변하여 종방향 항법 성능이 달라지기 때문이다. Fig. 4의 우측 그래프 보면 진행 방향이 Y축으로 바뀌는 14s에서 46s 사이에 종방향 표준 편차(점선)가 감소하며, 그 결과 위치 오차 경계가 감소하게 된다.

Fig. 4. Along-track position error bound results for nominal and GNSS failure condition (left) and ECEF standard deviation (right) of the GPS only case.

Fig. 5는 GPS와 Galileo를 함께 사용하는 상황에 대한 시뮬레이션 결과를 나타낸다. Fig. 5의 좌측 그래프를 보면, 다중 위성군의 위치 오차 경계는 단일 위성군의 경우와 같이 점진적으로 증가하고 단일 고장보다 이중 고장에서 더 큰 것을 확인할 수 있다. 하지만, 다중 위성군의 경우 가시 위성 수가 2배 정도로 늘어 위성 기하가 개선되면서 Fig. 4의 결과보다 위치 오차의 표준 편차가 감소한 것을 확인할 수 있다. 따라서, 오차 경계는 단일과 이중 고장 각각 5 m와 10 m로 수준으로 단일 위성군 결과에 비해 크게 감소한다. 또한 진행 방향에 따른 변화가 단일 위성군에 비해 작아진다.

Fig. 5. Along-track position error bound results for nominal and GNSS failure condition (left) and ECEF standard deviation (right) of the GPS+Galileo case.

4. CONCLUSIONS

본 논문에서는 위치 오차 경계 기반 geo-fencing function에 활용될 속도 불확실성을 고려한 정상상태 및 GNSS 고장상태에 대한 위험성 확률 평가 기반 종방향 위치 오차 경계를 산출 방법론을 제안하였다. 정확한 속도 추정을 위해 의사거리 측정치와 의사거리 변화율 측정치를 함께 사용하는 경우를 고려하였으며, 의사거리 변화율 측정치 고장의 영향을 고려한 위치 오차 경계 산출 방법론을 제안하였다. 이를 바탕으로 종방향 위치 오차 경계 산출 시뮬레이션을 수행하였다. 항법 시스템의 성능이 오차 경계에 미치는 영향을 분석하기 위해 고장 위성수와 위성군에 따른 결과를 비교하였다. 그 결과 종방향 위치 오차 경계는 기체의 자세에 크게 의존하는 것을 확인할 수 있었으며, 다중 위성군을 사용할 경우 단일 위성군의 40% 수준으로 줄일 수 있었다. 본 논문에서 제안된 고장상태에 대한 위험성 확률 평가 기반 위치 오차 경계 산출 방법론은 고장의 영향을 정량화할 수 있는 센서의 경우 해당 센서의 고장상태에 대한 안전성 요구조건을 만족하는 위치 오차 경계 도출에 동일하게 적용될 수 있다. 또한 임무 경로가 지정된 경우 비행기술오차를 함께 고려하여 총 시스템 오차에 대한 위치 오차 경계 산출로 확장이 가능하다. 따라서 향후 SAFEGAURD 뿐 아니라 안전성을 만족하는 UAS의 공역 설계, 충돌 회피 및 경로 설정 등에도 활용이 가능할 것으로 기대된다.

APPENDIX A

Joerger & Pervan (2013)에 따르면, i번째 고장 가설 하에 최악 상황 고장을 찾기 위해 식 (15)의 기울기를 식 (17)의 고장 벡터에 대한 식으로 표현한다. 식 (7, 10, 16)을 활용해 추정치 오차와 이노베이션 벡터의 기댓값을 나타내면 식 (A-1, A-2)와 같다.

\(E[\tilde{x}_k^+ ]=(I-K_k H_k ) F_{k-1} E[\tilde{x}_{k-1}^+ ]+K_k T_i f_{i,k}=K_k^" F_{k-1} E[\tilde{x}_{k-1}^+ ]+K_k T_i f_{i,k}\)                         (A-1)

\(E[\gamma_k  ]=-Η_k F_{k-1}  E[\tilde{x}_{k-1}^+  ]+T_i f_{i,k}\)                                                                     (A-2)

식 (A-1)에 따라 1부터 \(k\)번째 에폭까지 오차를 누적해 \(E[\tilde{x}_k^+ ]\)를 정리하면 식 (A-3, A-4)와 같다.

\(E[\tilde{x}_k^+]=[A_{1k}…A_{kk}] \left[\begin{array}{c}T_i f_{i,1}  \\⋮\\ T_i f_{i,k}  \end{array}\right]=A_{1:k}  f_{i,1:k}\)                                                           (A-3)

where    

\(A_{mk}=\left\{ \begin{array}{c}K^"_k K^"_{k-1}…K^"_{m+1} K_m, ~{\rm if}~ m<k \\ K_m, \end{array}\right.\)                                                         (A-4)

식 (A-3)을 (A-2)에 대입하여 \(E[\gamma_k ]\)를 구하면 식 (A-5)와 같다.

\(E[\gamma_k ]=[-H_k F_{k-1} A_{1:k-1} ~~ I] \left[\begin{array}{c}T_i f_{i,1} \\ ⋮ \\ T_i f_{i,k} \end{array}\right]=B_k  f_{i,1:k}\)                                                   (A-5)

식 (A-5)를 활용해 검정 통계량의 non-centrality를 구하면 식 (A-6, A-7, A-8)과 같다.

\(\lambda_k^2=\sum\limits_{j=1}^k f_{i,1:j}^T Τ_{i,j}^T B_j^T S_j^{-1}  B_j Τ_{i,j} f_{i,1:j}\)                                                                                                                

\(=f_{i,1:k}^T T_{i,k}^T [\bar{B}_1^T… \bar{B}_k^T ] \left[\begin{array}{ccc}S_1^{-1}& & \\ &⋱& \\ & &S_k^{-1}\end{array}\right] \left[\begin{array}{c}\bar{B}_1\\⋮\\\bar{B}_k \end{array}\right] Τ_{i,k} f_{i,1:k}\)                                                      

\(=f_{i,1:k}^T Τ_{i,k}^T \bar{B}_{1:k}^T S_{1:k}^{-1} \bar{B}_{1:k} T_{i,k} f_{i,1:k}\)                                                                                                       (A-6)

where    

\(T_{i,j}=\left[\begin{array}{ccc}T_i&0&0 \\ 0&⋱&0 \\ 0&0&T_i \end{array}\right]\in \mathbb{R}^{2Nj\times 2j}\)                                                                 (A-7)

\(\bar{B}_i=[(B_i ~~ 0_{2N\times 2N(k-i)}]_{2N\times 2N}\)                                                                     (A-8)

여기서 \(N\)은 가시 위성 수이다. 식 (A-3)과 (A-6)을 식 (15)에 대입하여 기울기를 고장 벡터에 대해 정리하면 식 (A-9)와 같다.

\(\alpha_{i,k}^2={f_{1:k}^T Τ_{i,k}^T A_k^T \tau_ε^T τ_ε A_k T_{i,k} f_{1:k}~~~~~~~~ \over f_{1:k}^T Τ_{i,k}^T \bar{B}_{1:k}^T S_{1:k}^{-1} \bar{B}_{1:k} T_{i,k} f_{1:k} }~~~~~~ \)                                                                 (A-9)

다음의 변수 변환을 정의하여 식 (A-9)를 정리하면,

\(\breve{f}_{1:k}=(Τ_{i,k}^T \bar{B}_{1:k}^T S_{1:k}^{-1} \bar{B}_{1:k} T_{i,k} )^{1/2} f_{1:k}\)                                                             (A-10)

\(M_z^{-1}=(Τ_{i,k}^T \bar{B}_{1:k}^T S_{1:k}^{-1} \bar{B}_{1:k} T_{i,k} )^{1/2}\)                                                                 (A-11)

\(M_x=\tau_\varepsilon A_k T_{i,k}\)                                                                                                 (A-12)

식 (A-9)는 최종적으로 식 (A-13)과 같은 내적꼴로 나타낼 수 있다.

\(\alpha_{1,k}^2={\breve{f}_{1:k}^T M_z^T M_x^T M_x M_z \breve{f}_{1:k}\over \breve{f}_{1:k,j}^T \breve{f}_{1:k,j} }=\left|M_x M_z\cdot {\breve{f}_{1:k} \over ‖\breve{f}_{1:k} ‖ } \right|^2\)                                                         (A-13)

따라서, 식 (A-13)은 \(M_x M_z\)와 \({\breve{f}_{1:k} \over ‖\breve{f}_{1:k} ‖ }\)가 평행할 때 기울기가 최대가 되며, 이때 최대 기울기(\(\alpha_{i,k,max}\))는 식 (A-14)와 같다.

\(\alpha_{i,k,max}=|M_x M_z |\)                                                                                         (A-14)

식 (A-14)는 한 고장 가설 하에서 최대 기울기이며, 고장 가설 중 가장 큰 기울기를 갖는 고장 가설이 최악 상황 고장 가설이 된다.