[Smoothing Algorithm] hatch filter

풍문으로만 들어봤던 Hatch Filter, 드디어 나왔다. 

 

GNSS 시스템을 다루다보면 다양한 필터를 만난다. 필자는 전공이 기계공학이 아니어서 제어에 약하다. 그러니 조금 더 꼼꼼하고 자세히 살펴보겠다. GNSS 데이터의 노이즈를 줄이기 위해 특화된 필터인 Hatch Filter에 대해 알아보자 

 

전리층을 통과하는 전파 신호의 지연으로 인해 코드 측정치, 반송파 위상 측정치 간의 차이가 벌어지는 현상인 code-carrier divergence이 발생한다. 전리층을 통과할 때 코드는 느려지고 반송파는 빨라져서 이 둘의 격차가 전리층 오차의 2배가 된다, 이러한 오차를 줄이는 것이 바로 Hatch Filter이다. 

 

 

---

Hatch Filter: 코드 - 반송파 평활화 

GPS 관측값에서 노이즈가 심한 코드(Code, C1)데이터를 매우 정밀한 위상(Carrier Phase, L1)데이터를 이용해 부드럽게 깎아주는(Smoothing)알고리즘 

 

해치 필터의 목적은 역시나 관측값의 노이즈를 줄여주는 것이다. 절대적인 위치는 코드가 잡고, 미세한 변화는 L1이 잡는 개념을 사용한다.

 

코드(C1)은 오차가 수 미터(m)로 크고 노이즈(Multi Path 등)가 심하지만 기준이 되는 절대적인 거리 정보는 정확하다.

반면 위상(L1)은 mm단위로 매우 정밀하지만, 처음에 미지정수(Ambiguity)라는 한계가 존재한다. 따라서 시작점의 절대적인 위치는 모르지만 직전 시간 대비 얼마나 움직였는지에 대한 변화량은 정밀하게 측정한다. 

 

해치 필터는 이 둘의 장점만 섞은 가중 평균(Weighted Average) 수식을 사용한다. 조금 더 쉽게 표현해보자면, 코드와 위상의 차이(평균 오차)를 실시간으로 계속 업데이트해서 위상에 더해주는 것과 같다. 

 

 

그럼 수식으로 알아보자 

• s: satellite(위성)
• k: k-th epoch(시점) 
• R: code measurement - 관측한 코드 의사거리 
• Φ: carrier measurement - 관측한 위상 데이터 
• R̂: smoothed code - 필터를 거친 뒤 부드러워진 최종 코드 결과값 

 

 

코드(R) 관측값은 1/n만큼 믿고, 이전 결과에 위상 변화량( Φ (s;k) - Φ (s;k-1))을 더한 추정값을 (n-1)/n만큼 크게 믿고 가중 평균을 내는 기본 식이다. 

 

다음은 해치 필터의 식을 수학적으로 이항하고 묶어서 맨 마지막 줄의 결론을 이끌어 내는 과정이다. 

 

3-4번째줄로 이어가는 과정을 치환했다고 생각하면 이해가 편하다. <R- Φ>는 1번째부터 k번째 시간까지 누적된 코드-위상의 최종 평균값을 의미한다. 새로운 식이 아니라, 현재까지의 평균이라는 기호로 표시한 것이다. 

 

그런데 이 알고리즘에는 치명적인 한계가 존재한다. 바로 위상 데이터에 사이클 슬립이 발생할 때마다 이 알고리즘은 초기화(n=1부터 다시 시작)해야한다. 

 

 

사실 개념 나름 이해할만하다. 코드와 반송파의 장점을 모아 코드로 초기 정확도를 잡고 추후 변동성을 반송파로 정밀하게 잡아나가는 것, 그런데 어려운 것은 이를 실제로 활용하고 다른 개념과 함께 사용하는 것 .. 

---

 

MATLAB을 통해 구현한 Hatch Filter + 이중 주파수 Multi Path 계산 

GPS L1/L2 이중주파수 멀티패스(Multipath) 계산과 Hatch Filter를 구현해보자. 

여기서 나오는 3.09라는 상수는 GPS 주파수에서 유도되는 고정 상수이다. 

 

•  GPS 신호 주파수
  • L1: f₁ = 1575.42 MHz
  • L2: f₂ = 1227.60 MHz

이중주파수를 이용해 코드(의사거리)와 반송파(위상)을 결합할 때 표준적으로 쓰이는 멀티패스 지표(MP1)의 공식은 다음과 같다. 

멀티패스(Multipath)란 ?

 

1. 물리적 현상으로서의 멀티패스 

위성 신호가 건물, 지면, 나무 등에 반사되어 수신기에 직접파(direct path)와 반사파(reflected path)가 섞여 들어오는 물리 현상으로 코드 추적(code tracking)과 반송파 추적(carrier tracking) 양쪽에 오차를 만듦 

 

2. 코드에서 다루는 MP1이라는 관측 지표(observable)

이 위성/이 시간대 신호 품질이 어떤가?를 나타내는 지표 

 

MP를 나타내는 이 식은 물리적 멀티패스만 뽑아낸 값이 아닙니다.

 

MP1 = P1 − (f₁²+f₂²)/(f₁²−f₂²)·L1 + 2f₂²/(f₁²−f₂²)·L2

 

이건 코드 의사거리와 반송파 위상의 차이에서, 기하학적 거리와 이온층 성분을 제거하고 남은 잔차(residual)이다. 이 잔차 안에는 사실 여러 가지가 섞여 들어간다. 

 

① 코드 멀티패스	반사파 효과 (코드 추적에 영향)
② 코드 수신기 잡음	열잡음, 양자화 오차 등 (thermal noise)
③ 반송파 위상의 정수 모호정수 바이어스	L1, L2 위상은 원래 미지의 정수 사이클(N)을 포함하는데, 이게 상수 바이어스로 남아있음 (시간에 따라 변하지 않으면 문제없지만, 슬립이 있으면 계단형으로 튐)
④ 코드-반송파 간 하드웨어 바이어스(DCB 등)	위성/수신기 장비 특성에 따른 차이
⑤ (반송파 자체의 멀티패스)	반송파도 미세하게 멀티패스 영향을 받지만 코드보다 훨씬 작음 →  파장 단위라서

즉 MP1 = 코드 멀티패스 + 잡음 + 일정한 바이어스 상수 를 의미한다. 

반송파(L1, L2)는 파장이 짧고 위상 추적이 정밀해서(mm~cm 수준) 상대적으로 잡음 없는 기준선(reference) 역할을 하고, 거기서 코드값을 빼면 코드에 있는 잡음+반사파 영향이 두드러지게 드러나는 것이다. 

 

MP1 = P1 − (f₁²+f₂²)/(f₁²−f₂²)·L1 + 2f₂²/(f₁²−f₂²)·L2

 

여기서 계수를 숫자로 계산해보면 

 

• f₁² = 2,481,948.18
• f₂² = 1,507,001.76
• f₁² − f₂² = 974,946.42
• 2f₂²/(f₁²−f₂²) = 3,014,003.52 / 974,946.42 ≈ 3.0913 ≈ 3.09
• (f₁²+f₂²)/(f₁²−f₂²) = 3,988,949.94 / 974,946.42 ≈ 4.0913 = 1 + 3.09

 

즉, 위 공식을 그대로 풀어쓰면:

MP1 = C1 − (1+3.09)·L1 + 3.09·L2 = C1 − L1 − 3.09(L1−L2)

 

아래 코드에서 3.09를 찾아볼 수 있는데, 이는 GPS L1/L2 반송파 주파수 비율에서 수학적으로 결정되는 상수인 것이다. 

 

 

% gLAB을 통해 읽어온 MEAS 파일 [C1C, L1C, C1P, L1P, C2P, L2P] 정보가 담김 
[~,~,~,~,~,~,~,dat_meas] = read_glab_v6(out_file);

for i = 1:N_prn
    prni = prns(i);
    ind_prn = dat_meas(:,4) == prni;
    sec_prn{i} = dat_meas(ind_prn,3);
    N_meas = length(sec_prn{i});

    % calculate MultiPath
    % MP_C1 = C1 - L1 - 3.09(L1-L2)
    C1 = dat_meas(ind_prn,8);
    L1 = dat_meas(ind_prn,11);
    L2 = dat_meas(ind_prn,13);
    MP_C1{i} = C1 - L1 - 3.09*(L1-L2);

    Ts = 60;
    C1s_2 = 0; L1p = 0;  %C1 smoothed
    C1s_3 = 0; L1fp = 0; %L1f previous

    MP_C1s_60_cur = zeros(N_meas,1);
    C1_DF_cur = zeros(N_meas,1);


    % Hatch Filter(C1 smoothing)
    for j = 1:N_meas
        % window size n
        n = min(j,Ts);

        L1f = L1(j) + 3.09*(L1(j) - L2(j));

        if j == 1
            C1s_2 = C1(j);
            C1s_3 = C1(j);
        else
            C1s_2 = C1(j)/n + ((n-1)/n) * (C1s_2 + (L1(j) - L1p));
            C1s_3 = C1(j)/n + ((n-1)/n) * (C1s_3 + (L1f  - L1fp));
        end

        % next Loop -> L1
        L1p = L1(j);
        L1fp = L1f;

        %Final Multipath
        MP_C1s_60_cur(j) = C1s_2- L1(j) -3.09*(L1(j)-L2(j));
        C1_DF_cur(j) = C1s_3 - L1f;

    end
    MP_C1s_60{i} = MP_C1s_60_cur;
    C1_DF{i} = C1_DF_cur;
end

 

 

변수	의미
MP_C1	평활화 없이 순간값으로 계산한 코드 멀티패스(표준 MP1)
MP_C1s_60	원시 L1로 60초 Hatch filter 평활화를 적용 후, 최종적으로 3.09 보정을 적용한 멀티패스
C1_DF	이온층 프리 반송파(L1f)로 60초 Hatch Filter를 평활화한 divergence-free 코드 - 반송파 결과

 

 

---

일단 아 Hatch Filter는 현재 시점까지 평활화되는 창(Window)이 있고 데이터가 쌓일수록 반송파로 추적되는 과거 이력의 비중이 커진다. 처음에는 코드 위주로 가다가 시간이 지날수록 반송파 위주로 전환되는 구조이다. 

 

이 예제에서 n = 60으로 두었지만 더 큰 값을 가질 수도 있다. 

 

일단 오늘은 여기까지 ! 
Multi Path에 대해서도 다시 정리한 뒤, 글을 작성해보겠다.