Modeling of Stochastic Process Noises for Kinematic GPS Positioning

GPS 이동측위를 위한 프로세스 잡음 모델링

Abstract

The Kalman filter has been widely used in the kinematic GPS positioning due to its flexibility and efficiency in computational points of view. At the same time, the relative positioning technique also provided the high precision positioning results by removing the systematic errors in the measurements significantly. However, the positioning quality may be degraded following to longer in baseline length. For this case, it is required that the remaining atmospheric effects, such as double-difference ionospheric delay and zenith wet delay, should be properly modeled by examining the characteristics of the stochastic processes. In general, atmospheric effects are estimated with the assumption of random walk, or the first-order Gauss-Markov stochastic process, which requires the precise modeling on the corresponding process noises. Therefore, we determined and provided the parameters for modelling the process noises for atmospheric effects. The auto-correlation functions are empirically determined at first, and then the parameters are extracted from the empirical auto-correlation function. In fact, the test results can be either applied directly, or used as guidance values for the modeling of process noises in the kinematic GPS positioning.

알고리즘의 유연성 및 효율성으로 인해 GPS 이동측위 시 칼만필터가 주로 사용되어 왔으며 동시에 다양한 계통오차의 제거가 가능한 상대측위 기법이 널리 사용되어 왔다. 하지만 기선의 길이가 길어지게 되면 상대측위 기법을 사용하더라도 대기효과를 충분히 제거하기 어렵기 때문에 이 경우 제거되지 않고 남아 있는 대기효과를 상태벡터에 추가하여 추정을 하기도 한다. 칼만필터를 이용하는 경우 일반적으로 대기효과는 랜덤워크 혹은 일차가우스-마르코프 프로세스로 모델링하게 되는데 이때 프로세스 잡음에 대한 정확한 모델링이 필수적이다. 본 연구에서는 대기효과에 해당되는 프로세스 잡음 모델링을 위해 필요한 매개변수를 결정하였다. 이를 위해 이중차분 전리층 지연값과 천정방향 습윤지연값을 이용하여 실험적 자기상관함수를 계산하였으며 이를 통해 프로세스 잡음 모델링에 필요한 매개변수를 계산하였다. 결정된 매개변수값들은 유사한 대기환경에서 취득된 데이터에 대한 프로세스 잡음 모델링 시 직접 사용될 수 있으며 유사한 대기환경이 아닌 경우일 지라도 초기 근사값으로 활용될 수 있을 것이다.

1. 서 론

GPS를 이용하여 정확한 위치결정을 위한 대표적인 방법은 상대측위기법으로 기준국과 이동국에서 동시에 수신된 데이터를 서로 빼줌(차분)으로써 대류권 지연과 전리층 지연 등을 포함한 계통오차의 소거가 가능하다(Hofmann-Wellenhof et al., 2001). 하지만 수신국과 이동국 사이의 거리 즉, 기선의 길이가 증가하게 되면 차분법을 사용하더라도 대기효과의 충분한 제거가 사실상 불가능하다. 따라서 일반적으로 중·장기선에 해당되는 데이터의 처리를 위해서는 남아있는 대기효과를 미지수(unknown)로 모델링하여 추정을 한다(Goad and Yang, 1997; Yang et al., 2000; Schüler, 2001). 중 · 장기선 기선처리는 기준점 측량과 같이 특정 지점의 정밀한 위치를 결정하는 데 사용되거나 이동하는 항공기 등의 위치를 결정하는 항법에 사용되기도 한다. 전자의 경우를 GPS 측위기법 중정지측위기법, 후자는 이동측위기법이라 하며 측지 · 측량 분야에는 분야 특성상 정지측위기법이 주로 사용된다. 하지만 항공중력측량과 같이 이동측위기법이 필요한 경우도 있으며 이 경우 항공기의 위치를 정확하게 결정하는 것이 매우 중요하다(Jekeli and Garcia, 1997). 이동측위를 위한 추정 방법은 다양하게 존재할 수 있으나 칼만필터(Kalman filter)가 실시간 데이터 처리 능력 및 처리 알고리즘의 효율성 등을 고려했을 때 가장 적합한 것으로 알려져 있다. 칼만필터는 크게 프로세스 방정식(process equation)과 관측방정식(observation equation)으로 구성되어 있다(Gelb, 1974; Brown and Hwang, 1997; Simon, 2006). 프로세스 방정식은 상태벡터의 시간에 따른 변화특성을 설명하는 부분으로 상태벡터를 어떤 종류의 프로세스로 정의하느냐와 해당 프로세스 잡음은 어떻게 설정할 것인가가 필수 요소이다. 대기효과에 해당되는 프로세스 모델도 다양하게 존재하나 중·장기선에 적합한 대기효과 추정 모델로는 랜덤워크(random walk)나 일차 가우스-마르코프(first-order Gauss-Markov) 모델이 주로 제안되었다(Schüler, 2001; Yang et al., 2000). 이유는 이 두 모델이 구현이라는 측면에서 비교적 간단하며 대기효과의 시간에 따른 변화를 충분히 표현할 수 있는 것으로 알려져 있기 때문이다. 하지만 두 가지 프로세스에 해당되는 잡음은 수신 환경 등에 영향을 받기 때문에 이에 대한 정확한 모델링에는 여전히 한계가 있다. 이러한 이유로 실험적인 방법을 통한 모델링을 위한 연구가 다양하게 진행되어 왔다(Gelb, 1974; Jekeli, 2001; Schüler, 2001). 본 연구에서는 이중차분(DD: Double-Difference) 전리층 지연과 천정방향 습윤지연(ZWD: Zenith Wet Delay)에 대한 프로세스 잡음 분석을 통해 유사한 환경에서 수신된 이동측위용 GPS 데이터 처리에 유용한 모델링 값을 제공하고자 하였다. 이를 위해 DD 전리층 지연과 ZWD를 이용하여 실험적(empirical) 자기상관함수(auto-correlation function)를 결정하였으며 이로부터 프로세스 잡음 모델링에 필요한 매개변수(parameter)를 계산하였다.

 

2. 이론적 고찰 및 연구방법

프로세스 잡음 특성은 이미 알고 있는 잡음 특성 정보를 활용하거나 초기 결과값 형태로 추정된 상태벡터값을 이용하여 결정할 수 있다. 하지만 잡음 특성에 대한 정보를 미리 알 수 있는 경우는 제한적이며 따라서 일반적으로는 추정된 상태벡터 값으로부터 잡음 특성에 대한 분석을 통해 잡음 특성을 재정의 하게 된다. 잡음 특성은 이미 알려져 있는 자기상관함수로부터 정의할 수 있으며 상태벡터에 대한 자기상관함수 정보가 없는 경우에는 실험적인 방법을 통해 자기상관함수를 결정하게 된다. 즉, 잡음 특성을 설명하는 매개변수는 실험적으로 결정된 자기상관함수로부터 추출할 수 있다.

2.1 프로세스 방정식과 잡음

칼만필터는 시간에 따라 변하는 양을 추정하는 기법 중에 하나이며 Eq. (1)에서와 같이 상태벡터(state vector)의 시간에 따른 변화 특성을 정의하는 프로세스 방정식과 상태벡터와 관측값 사이의 관계를 정의하는 관측방정식으로 구성되어 있다(Gelb, 1974; Brown and Hwang, 1997).

where x: state vector, Φ: state transition matrix, w: process noise, z: measurements, H: design matrix, v: measurement noise, Q : covariance matrix for process noise, R : covariance matrix for measurement noise.

GPS 관측값에 포함되어 있는 대기효과는 크게 대류권 지연과 전리층 지연으로 나눌 수 있으며 각각에 해당되는 상태 벡터의 시간에 따른 변화는 랜덤워크 또는 일차 가우스-마르코프 프로세스로 모델링 할 수 있다.

랜덤워크를 나타내는 이산 프로세스(discrete process)는 Eq. (2)와 같으며 해당 프로세스 잡음에 대한 분산 qk는 스펙트럼 강도(spectral intensity)라 불리는 q와 시간간격(time interval)인 Δt의 곱으로 표현되며 Eq. (3)과 같이 나타낼 수 있다(Gelb , 1974).

Eq. (2)에서 알 수 있듯이 이전 epoch에서의 상태벡터에 백색잡음(white noise)이 더 해진 형태로 다음 epoch에서의 상태벡터가 결정된다. 일차 가우스-마르코프 프로세스의 경우 Eq. (4)에서 볼 수 있듯이 다음 epoch에서의 상태벡터 계산을 위해서는 상관시간(correlation time)의 역수인 β에 대한 정보가 필요하다. 뿐만 아니라 프로세스 잡음에 대한 분산의 계산하기 위해서는 상태(state)에 대한 분산 정보도 필요하다.

즉, 정확한 프로세스 잡음 모델링을 위해서는 스펙트럼 강도(q), 상관시간(1/β), 상태에 대한 분산 값들을 결정하여야 한다

2.2 프로세스 잡음 모델링 파라미터 결정

프로세스 잡음 모델링을 위한 파라미터는 자기상관함수에 대한 정보가 있는 경우 이를 이용하여 구할 수 있으나 자기상관함수에 대한 정보가 없는 경우가 일반적이기 때문에 이 경우 실제 데이터를 이용하여 실험적 자기상관함수를 결정하게 된다. Eqs. (6) and (7)는 실험적 자기상관함수를 계산하기 위한 식을 나타낸 것이다(Gelb, 1974).

where C : empirical auto-correlation function for x, m: mean of x, N: total number of x, Δt: time interval.

따라서 초기 결과로 획득된 대기효과에 의한 지연량인 DD 전리층 지연값과 ZWD에 대한 실험적 자기상관함수를 계산하여야 한다. 계산된 자기상관함수로부터 스펙트럼 강도, q는 Eq. (8)을 통해 계산할 수 있다(Schüler, 2001).

상관시간은 자기상관함수의 값이 분산값의 1/e에 해당되는 시간으로 정의하고 있기 때문에 실험적 자기상관함수가 결정되기만 하면 쉽게 계산이 가능하다.

 

3. 자료처리 및 결과 분석

전리층 지연 및 대류권 지연값에 대한 실험적 자기상관함수를 얻기 위해 DGIST(2012)가 제시한 방법을 이용하여 GPS 데이터를 처리한 후 이동측위 결과와 DD 전리층 지연 및 ZWD에 대한 초기값을 추출하였다. 2007년 1월 3일 0시부터 4시 18분 까지(UTC 기준) 김포(GIMP)에 설치된 기준국과 항공기(AIRP)에 탑재된 GPS 수신기를 통해 수신된 이중주파수 데이터를 이용하였으며 상대측위기법을 통해 DD 전리층 지연과 ZWD를 얻을 수 있었다. 데이터 취득 시 Kp-index값은 대략 3 정도였으며 따라서 전리층 지연의 변화는 심하지 않은 경우에 해당된다고 볼 수 있다. 왜냐하면 미국의 NOAA(National Oceanic and Atmospheric Adminstration)는 Kp-index가 4 이상일 때 지자기 교란이 있다고 판단하며 이때 경보를 발령한다. 참고로 2009년 한 해 동안 Kp-indx를 보면 3 이하인 경우가 98%를 차지하여 대부분의 경우 지자기 교란 정도가 심하지 않다는 것을 알 수 있다. Fig. 1은 GIMP 기준국의 위치와 GIMP를 출발하여 부산 근처까지 왕복한 항공기의 궤적을 표시한 것이다. 이때 이륙 및 착륙을 제외한 나머지 구간에서는 대략 3km의 고도를 유지하였다.

Fig. 1.Location of reference station and trajectories of aircraft

기준국은 GIMP 하나만 사용하였기 때문에 단일기선 처리에 해당되며 Fig. 2는 GIMP-AIRP 기선에 대한 DD 전리층 지연값을 시간에 따라 도시한 것이다. Fig. 2에서 보이는 연속 데이터는 각 위성쌍(satellite pair)에 해당되는 DD 전리층 지연값을 의미한다. 따라서 데이터가 단절된 epoch에서는 위성쌍의 조합이 달라졌다는 것을 의미한다. DD 전리층 지연값의 시간에 따른 변화가 상대적으로 다르게 나타나는 데 이는 위성의 기하학적 배치와 고각이 서로 다르기 때문이다. 값의 범위는 대략 -50cm ∼ +50cm로 시간에 따른 변화량 자체는 크지 않은 것으로 나타났다.

Fig. 2.DD ionospheric delays (GIMP-AIRP baseline)

Fig. 2의 데이터는 모든 위성쌍에 해당되는 DD 전리층 지연값이므로 각각의 위성쌍에 해당되는 데이터로부터 실험적 자기상관함수를 결정하여야 한다. 이 때 연속 데이터 구간이 짧을 경우 신뢰할 수 있는 결과를 얻기 어렵기 때문에 연속 데이터 구간이 1시간보다 큰 데이터구간만을 선택하였다. Fig. 3는 하나의 위성쌍에 해당되는 DD 전리층 지연값을 예제로 도시한 것이다. 이 경우 연속 데이터 구간이 대략 2시간 이상이며 대략 -25cm ∼ -5cm 사이의 값을 보이고 있다.

Fig. 3.Example of DD ionospheric delays

Fig. 3에 해당되는 실험적 자기상관함수는 Eqs. (5) and (6)를 통해 계산할 수 있으며 그 결과는 Fig. 4와 같다. 이 경우 자기상관함수로부터 분산은 자기상관함수의 첫 번째 데이터를 의미하기 때문에 2.0×10-3m2이며 상관시간(1/β)은 약 20분(1186초)으로 계산되었다.

Fig. 4.Empirically determined auto-correlation function for the data shown in Fig. 3

위와 같이 1시간 이상의 연속 데이터 구간을 갖는 DD 전리층 지연값을 갖는 경우는 총 11개로 조사되었으며 각각에 대해 실험적 자기상관함수를 마찬가지로 계산하여 도시하면 Fig. 5와 같다. Fig. 5에서 나타난 경향을 보면 두 데이터는 자기상관이 급격히 떨어지는 반면 나머지 데이터들은 비교적 완만한 하강을 보이고 있는 특징이 있다.

Fig. 5.Empirically determined auto-covariance functions for DD ionospheric delays

Fig. 5에 나타난 실험적 자기상관함수로부터 2절에서 제시한 방법으로 스펙트럼 강도, 분산, 상관시간을 계산할 수 있으며 그 결과를 요약하면 Table 1과 같다.

Table 1.Spectral intensities, variances and correlation times derived from the DD ionospheric delay

11개의 연속 데이터 구간을 통해 계산된 결과에 의하면 스펙트럼 강도는 7.7×10-6 m2/s, 분산은 1.8×10-3 m2, 상관시간은 약 14분(830초)의 평균값을 보였다. 앞에서 언급하였듯이 데이터 취득 시 전리층 교란 정도는 심하지 않은 경우이므로 유사한 전리층 교란 상황에서의 GPS 이동측위 시 프로세스 모델링 값으로 직접 사용하여도 무방할 것으로 판단된다. 뿐만 아니라 전리층 교란이 상대적으로 큰 경우일 지라도 본 실험 결과의 값을 참고하여 DD 전리층 지연에 대한 대략적인 스펙트럼 강도, 분산, 상관시간을 정의하는 데 도움을 줄 수 있을 것이다.

천정방향 대류권 지연은 건조지연과 습윤지연으로 구성되어 있으며 건조지연량은 비교적 정확한 예측이 가능하기 때문에 모델링된 값을 사용하여도 오차가 거의 없는 것으로 알려져 있다(Hofmann-Wellenhof et al., 2001). 이에 반해 ZWD은 수증기량의 시간에 따른 변화로 인해 그 크기는 비록 작은 편이지만 정확한 예측이 어려운 것으로 알려져 있다. ZWD를 보다 구체적으로 나타내면 Eq. (9)과 같으며 ZWD는 ZWD 모델값과 잔차의 합으로 나타낼 수 있다.

where ZWDmodel: zenith wet delay value calculated from the model, ∆ZWD: zenith wet delay residual.

따라서 데이터 처리를 통해 직접 ZWD를 계산하거나 습윤지연 모델값을 기준으로 한 잔차(∆ZWD: ZWD residual)를 추정하는 방법을 주로 사용하게 된다. 두 경우 프로세스의 특성이 다르게 나타날 수 있기 때문에 각각의 경우에 대한 분석을 수행하였다. 먼저 ZWD 잔차와 이 잔차에 ZWD 모델값을 더한 전체 ZWD를 계산하였으며 Fig. 6은 그 결과를 도시한 것이다. 따라서 Fig. 6에 나타난 값들은 항공기 위치에서의 ZWD를 나타낸 것이다.

Fig. 6.Variations of ZWD and ZWD residuals with respect to time

Fig. 6에서 보는 바와 같이 값의 범위는 5cm에서 30cm 사이로 ZWD와 ZWD 잔차 사이에는 편이(bias) 형태의 차이를 보이는 것으로 나타났다. 이 두 경우에 대해 실험적 자기상관함수를 DD 전리층 지연 경우와 마찬가지로 계산할 수 있으 며 그 결과는 Fig. 7과 같다. ZWD와 잔차에 대한 분산은 각각 2.2×10-3m2, 1.5×10-3m2로 나타났으며 상관시간은 약 36분(2180초), 49분(2962초)으로 각각 계산되었다.

Fig. 7.Empirically determined auto-correlation functions for ZWD and ZWD residuals

Table 2는 실험적 자기상관함수로부터 계산된 ZWD 및 잔차에 대한 스펙트럼 강도, 분산, 상관시간을 보여주고 있다. 스펙트럼 강도와 분산은 각각 10-7, 10-3 수준으로 계산되었으며 대류권 지연의 경우 DD 전리층 지연의 경우에 비해 상관시간이 길게 나타났다. 이로부터 대류권 지연의 경우가 시간에 따른 예측이 상대적으로 용이하다는 것을 알 수 있다.

Table 2.Spectral intensities, variances and correlation times derived from the ZWD and ZWD residuals

Table 1과 Table 2는 상대측위를 이용한 이동측위 결과로부터 추출된 DD 전리층 지연값과 ZWD를 이용하여 계산된 값이므로 모든 이동측위에 직접 적용하는 데에는 한계가 있을 수 있으나 유사한 대기환경으로 판단되는 경우 제시된 결과값을 직접 사용하여도 무방할 것으로 판단된다. 뿐만 아니라 유사한 대기환경이라고 판단하기 어려운 경우일지라도 본 연구에서 제시된 값을 근사치로 활용하여 최적의 값들을 재정의 하는데 이용될 수 있을 것으로 예상된다.

 

4. 결론 및 토의

본 연구에서는 칼만필터를 이용한 이동측위 시 대기효과인 DD 전리층 지연과 ZWD의 추정을 위한 모델링에 필요한 스펙트럼 강도, 분산 그리고 상관시간을 결정한 후 활용성을 제시하고자 하였다. 이를 위해 DD 전리층 지연과 ZWD를 이용하여 실험적 자기상관함수를 계산하였으며 이로부터 프로세스 잡음 모델링에 필요한 스펙트럼 강도, 분산, 상관시간을 결정하였다. 본 연구를 통해 얻은 결론을 요약하면 다음과 같다.

첫째, DD 전리층 지연값은 위성조합에 따라 다양한 형태의 변화를 보일 수 있으며 연속데이터 구간의 길이도 달라지게 된다. 따라서 1시간 이상 연속데이터 구간만을 선택하여 각 데이터 구간에 대해 스펙트럼 강도, 분산, 상관시간을 결정하였으며 이로부터 스펙트럼 강도는 7.7×10-6m2/s, 분산은 1.8×10-3m2, 상관시간은 14분(830초)의 값을 얻었다. 이 값들은 전리층 교란이 비교적 적은 대기환경(Kp-index≤3)하에서 취득된 데이터를 기반으로 계산되었기 때문에 유사한 대기환경에서 취득된 데이터를 처리하는 경우 본 연구를 통해 제시된 값을 DD 전리층 지연값의 추정에 직접 사용하여도 무방할 것으로 판단된다. 또한 전리층 교란이 심한 대기환경하에서는 적절한 근사값으로의 역할을 할 수 있을 것으로 예상된다.

둘째, 대류권 지연은 추정하는 형태에 따라 ZWD 혹은 ZWD 잔차 형태로 나눌 수 있으며 각각의 경우에 대해 스펙트럼 강도, 분산, 상관시간을 결정하였다. ZWD 형태로 자기상관함수를 계산한 경우 스펙트럼 강도는 6.7×10-7m2/s, 분산은 2.2×10-3m2, 상관시간은 36분(2180초) 나타났으며 ZWD 잔차의 형태로 계산 한 경우 각각 1.9×10-7m2/s, 1.5×10-3m2, 49분(2962초)으로 나타났다. 또한 상관시간은 대기환경에 따라 다소 차이를 보일 수 있지만 본 연구에서 제시한 값이 유의미한 수준의 근사값으로 활용될 수 있을 것으로 예상된다.