복붙노트

[PYTHON] 위치 데이터에 파이썬에서 칼만 필터를 사용하는 방법?

PYTHON

위치 데이터에 파이썬에서 칼만 필터를 사용하는 방법?

[편집하다] @Claudio의 대답은 이상 치를 필터링하는 방법에 대한 좋은 팁을 제공합니다. 그래도 내 데이터에 칼만 필터를 사용하고 싶습니다. 그래서 아래의 예제 데이터를 변경하여 너무 극한이 아닌 미세한 변이 노이즈를 발생 시켰습니다 (많은 것을 볼 수 있습니다). 다른 사람들이 내 데이터에 PyKalman을 사용하는 방법에 대한 지침을 줄 수 있다면 좋을 것입니다. [/편집하다]

로봇 프로젝트의 경우 카메라로 공중파 연을 추적하려고합니다. Python으로 프로그래밍 중이며 아래에 시끄러운 위치 결과를 붙여 넣었습니다 (모든 항목에도 datetime 객체가 포함되어 있지만 명확성을 위해 제외 시켰습니다).

[           # X     Y 
    {'loc': (399, 293)},
    {'loc': (403, 299)},
    {'loc': (409, 308)},
    {'loc': (416, 315)},
    {'loc': (418, 318)},
    {'loc': (420, 323)},
    {'loc': (429, 326)},  # <== Noise in X
    {'loc': (423, 328)},
    {'loc': (429, 334)},
    {'loc': (431, 337)},
    {'loc': (433, 342)},
    {'loc': (434, 352)},  # <== Noise in Y
    {'loc': (434, 349)},
    {'loc': (433, 350)},
    {'loc': (431, 350)},
    {'loc': (430, 349)},
    {'loc': (428, 347)},
    {'loc': (427, 345)},
    {'loc': (425, 341)},
    {'loc': (429, 338)},  # <== Noise in X
    {'loc': (431, 328)},  # <== Noise in X
    {'loc': (410, 313)},
    {'loc': (406, 306)},
    {'loc': (402, 299)},
    {'loc': (397, 291)},
    {'loc': (391, 294)},  # <== Noise in Y
    {'loc': (376, 270)},
    {'loc': (372, 272)},
    {'loc': (351, 248)},
    {'loc': (336, 244)},
    {'loc': (327, 236)},
    {'loc': (307, 220)}
]

처음에는 특이 값을 계산 한 다음 실시간으로 데이터에서 간단히 제거했습니다. 그런 다음 필자는 칼만 필터에 대해 읽고 어떻게 잡음이 많은 데이터를 원활하게 만들 것인지 구체적으로 설명했습니다. 그래서 어떤 검색을 한 후에 나는 이것을 위해 완벽 해 보이는 PyKalman 라이브러리를 발견했습니다. 필자는 전체 칼만 필터 용어에서 다소 잃어 버렸기 때문에 위키와 칼만 필터의 다른 페이지를 읽었습니다. 나는 칼만 필터에 대한 일반적인 아이디어를 얻었지만, 필자는 그것을 코드에 어떻게 적용해야하는지에 대해서는 정말로 분실했다.

PyKalman 문서에서 다음 예제를 발견했습니다.

>>> from pykalman import KalmanFilter
>>> import numpy as np
>>> kf = KalmanFilter(transition_matrices = [[1, 1], [0, 1]], observation_matrices = [[0.1, 0.5], [-0.3, 0.0]])
>>> measurements = np.asarray([[1,0], [0,0], [0,1]])  # 3 observations
>>> kf = kf.em(measurements, n_iter=5)
>>> (filtered_state_means, filtered_state_covariances) = kf.filter(measurements)
>>> (smoothed_state_means, smoothed_state_covariances) = kf.smooth(measurements)

나는 단순히 다음과 같이 내 관찰에 대한 관찰을 대체했다.

from pykalman import KalmanFilter
import numpy as np
kf = KalmanFilter(transition_matrices = [[1, 1], [0, 1]], observation_matrices = [[0.1, 0.5], [-0.3, 0.0]])
measurements = np.asarray([(399,293),(403,299),(409,308),(416,315),(418,318),(420,323),(429,326),(423,328),(429,334),(431,337),(433,342),(434,352),(434,349),(433,350),(431,350),(430,349),(428,347),(427,345),(425,341),(429,338),(431,328),(410,313),(406,306),(402,299),(397,291),(391,294),(376,270),(372,272),(351,248),(336,244),(327,236),(307,220)])
kf = kf.em(measurements, n_iter=5)
(filtered_state_means, filtered_state_covariances) = kf.filter(measurements)
(smoothed_state_means, smoothed_state_covariances) = kf.smooth(measurements)

그러나 그것은 저에게 의미심장 한 자료를주지 않는다. 예를 들어 smoothed_state_means는 다음과 같습니다.

>>> smoothed_state_means
array([[-235.47463353,   36.95271449],
       [-354.8712597 ,   27.70011485],
       [-402.19985301,   21.75847069],
       [-423.24073418,   17.54604304],
       [-433.96622233,   14.36072376],
       [-443.05275258,   11.94368163],
       [-446.89521434,    9.97960296],
       [-456.19359012,    8.54765215],
       [-465.79317394,    7.6133633 ],
       [-474.84869079,    7.10419182],
       [-487.66174033,    7.1211321 ],
       [-504.6528746 ,    7.81715451],
       [-506.76051587,    8.68135952],
       [-510.13247696,    9.7280697 ],
       [-512.39637431,   10.9610031 ],
       [-511.94189431,   12.32378146],
       [-509.32990832,   13.77980587],
       [-504.39389762,   15.29418648],
       [-495.15439769,   16.762472  ],
       [-480.31085928,   18.02633612],
       [-456.80082586,   18.80355017],
       [-437.35977492,   19.24869224],
       [-420.7706184 ,   19.52147918],
       [-405.59500937,   19.70357845],
       [-392.62770281,   19.8936389 ],
       [-388.8656724 ,   20.44525168],
       [-361.95411607,   20.57651509],
       [-352.32671579,   20.84174084],
       [-327.46028214,   20.77224385],
       [-319.75994982,   20.9443245 ],
       [-306.69948771,   21.24618955],
       [-287.03222693,   21.43135098]])

나보다 밝은 영혼이 나에게 올바른 방향으로 어떤 힌트 나 예를 줄 수 있을까요? 모든 팁을 환영합니다!

해결법

  1. ==============================

    1.TL, DR, 하단의 코드와 그림을 참조하십시오.

    TL, DR, 하단의 코드와 그림을 참조하십시오.

    저는 칼만 필터가 여러분의 어플리케이션에서 꽤 잘 작동한다고 생각합니다.하지만 카이트의 역학 / 물리학에 대해 좀 더 생각할 필요가 있습니다.

    이 웹 페이지를 읽는 것이 좋습니다. 나는 저자와 관련이 없거나 지식이 없지만 칼만 필터에 머리를 대고 하루를 보냈다.이 페이지는 정말로 나를 위해 그것을 만들었다.

    간단히; (즉, 상태와 입력을 알고 미래 상태를 예측할 수있는) 동적 인 시스템에 대해 시스템에 대해 알고있는 것을 결합하여 실제 상태로 추정하는 최적의 방법을 제공합니다. 영리한 비트 (설명하는 페이지에서 볼 수있는 모든 행렬 대수에 의해 처리됩니다)는 두 가지 정보를 최적으로 결합하는 방법입니다.

    당신은 당신이 각각의 (공분산 행렬 R과 Q를 통해) 각각의 상태를 어떻게 지정하고, 칼만 게인 (Kalman Gain)은 모델을 얼마나 믿어야하는지 (예 : 당신의 현재 상태 추정치)와 얼마나 많이 당신의 측정을 믿으십시오.

    더 이상 걱정하지 않고 연의 간단한 모델을 만들어 봅시다. 제가 아래에 제안하는 것은 매우 간단한 모델입니다. 카이트의 역 동성에 대해 더 많이 알고 있으므로 더 나은 것을 만들 수 있습니다.

    연을 입자로 다루겠습니다 (분명히 단순화, 실제 연은 확장 된 몸체이므로 3 차원으로 방향이 있습니다). 편의를 위해 상태 벡터에 쓸 수있는 4 가지 상태가 있습니다.

    x = [x, x_dot, y, y_dot],

    여기서 x와 y는 위치이고 _dot는 각 방향의 속도입니다. 귀하의 질문에서 나는 우리가 측정 벡터에 쓸 수있는 두 가지 (잠재적으로 시끄러운) 측정이 있다고 가정합니다 :

    z = [x, y],

    측정 행렬 (여기에서 논의 된 H와 pykalman 라이브러리의 관측 행렬)을 기록 할 수 있습니다.

    z = Hx => H = [[1, 0, 0, 0], [0, 0, 1, 0]]

    그런 다음 시스템 역학을 설명해야합니다. 여기서 나는 외부 힘이 작용하지 않으며 카이트의 움직임에 감쇠가 없다고 가정 할 것입니다. (지식이 많을수록 더 잘 할 수있을 것입니다. 이것은 효과적으로 외력과 댐핑을 알려지지 않은 / 조절되지 않은 교란으로 간주합니다).

    이 경우, 이전 샘플 "k-1"의 상태의 함수로서 현재 샘플 "k"에서 각 상태에 대한 동역학은 다음과 같이 주어진다 :

    x (k) = x (k-1) + dt * x_dot (k-1)

    x_dot (k) = x_dot (k-1)

    y (k) = y (k-1) + dt * y_dot (k-1)

    y_dot (k) = y_dot (k-1)

    여기서 "dt"는 시간 단계입니다. 현재 위치와 속도를 기반으로 (x, y) 위치가 업데이트되고 속도는 변경되지 않은 것으로 가정합니다. 단위가 주어지지 않는다면 위의 방정식, 즉 position_units / sample_interval 단위로 "dt"를 생략 할 수있는 속도 단위라고 말할 수 있습니다 (측정 된 샘플이 일정한 간격으로 있다고 가정합니다). 이 4 개의 방정식을 역학 행렬로 요약 할 수 있습니다 (여기서는 F, pykalman 라이브러리의 transition_matrices).

    [0,1,0,0], [0, 0, 1, 1], [0, 0], F (k-1) = F = [[1, 1, 0, 0] , 0, 1]].

    우리는 이제 Python에서 Kalman 필터를 사용할 수 있습니다. 귀하의 코드에서 수정 된 항목 :

    from pykalman import KalmanFilter
    import numpy as np
    import matplotlib.pyplot as plt
    import time
    
    measurements = np.asarray([(399,293),(403,299),(409,308),(416,315),(418,318),(420,323),(429,326),(423,328),(429,334),(431,337),(433,342),(434,352),(434,349),(433,350),(431,350),(430,349),(428,347),(427,345),(425,341),(429,338),(431,328),(410,313),(406,306),(402,299),(397,291),(391,294),(376,270),(372,272),(351,248),(336,244),(327,236),(307,220)])
    
    initial_state_mean = [measurements[0, 0],
                          0,
                          measurements[0, 1],
                          0]
    
    transition_matrix = [[1, 1, 0, 0],
                         [0, 1, 0, 0],
                         [0, 0, 1, 1],
                         [0, 0, 0, 1]]
    
    observation_matrix = [[1, 0, 0, 0],
                          [0, 0, 1, 0]]
    
    kf1 = KalmanFilter(transition_matrices = transition_matrix,
                      observation_matrices = observation_matrix,
                      initial_state_mean = initial_state_mean)
    
    kf1 = kf1.em(measurements, n_iter=5)
    (smoothed_state_means, smoothed_state_covariances) = kf1.smooth(measurements)
    
    plt.figure(1)
    times = range(measurements.shape[0])
    plt.plot(times, measurements[:, 0], 'bo',
             times, measurements[:, 1], 'ro',
             times, smoothed_state_means[:, 0], 'b--',
             times, smoothed_state_means[:, 2], 'r--',)
    plt.show()
    

    다음은 노이즈를 거부하는 합리적인 작업을 보여주는 다음과 같은 결과입니다 (파란색은 x 위치, 빨간색은 y 위치, x 축은 샘플 번호 임).

    위의 그림을보고 너무 울퉁불퉁하다고 생각해 봅시다. 어떻게 고칠 수 있니? 위에서 설명한 것처럼 칼만 필터는 두 가지 정보에 작용합니다.

    위의 모델에서 얻은 다이내믹은 매우 간단합니다. 문자 그대로 말하면, 위치는 현재 속도 (명백하고 물리적으로 합리적인 방식으로)에 의해 업데이트되고 속도는 일정하게 유지됩니다 (이것은 분명히 물리적으로 사실이 아니지만 속도가 천천히 변화해야한다는 직관을 포착합니다).

    추정 된 상태가 더욱 부드럽게되어야한다고 생각하면,이를 달성하는 한 가지 방법은 우리의 역학보다 측정에 대한 신뢰도가 낮다는 것입니다. 즉, 우리는 state_covariance에 비해 관측치의 편차가 더 큽니다.

    위 코드의 끝에서 시작하여 관측 공분산을 이전에 추정 된 값의 10 배로 고정하고, 관찰 공분산의 재평가를 피하기 위해 em_vars를 그림과 같이 설정해야합니다 (여기를 참조하십시오)

    kf2 = KalmanFilter(transition_matrices = transition_matrix,
                      observation_matrices = observation_matrix,
                      initial_state_mean = initial_state_mean,
                      observation_covariance = 10*kf1.observation_covariance,
                      em_vars=['transition_covariance', 'initial_state_covariance'])
    
    kf2 = kf2.em(measurements, n_iter=5)
    (smoothed_state_means, smoothed_state_covariances)  = kf2.smooth(measurements)
    
    plt.figure(2)
    times = range(measurements.shape[0])
    plt.plot(times, measurements[:, 0], 'bo',
             times, measurements[:, 1], 'ro',
             times, smoothed_state_means[:, 0], 'b--',
             times, smoothed_state_means[:, 2], 'r--',)
    plt.show()
    

    아래 그림을 생성합니다 (점선으로 측정, 점선으로 상태 추정). 그 차이는 다소 미묘하지만, 더 부드럽다는 것을 알 수 있기를 바랍니다.

    마지막으로,이 맞는 필터를 온라인으로 사용하려면 filter_update 메소드를 사용하면됩니다. 매끄러운 방법은 배치 측정에만 적용 할 수 있기 때문에이 방법은 매끄러운 방법이 아닌 필터 방법을 사용합니다. 여기 더 :

    time_before = time.time()
    n_real_time = 3
    
    kf3 = KalmanFilter(transition_matrices = transition_matrix,
                      observation_matrices = observation_matrix,
                      initial_state_mean = initial_state_mean,
                      observation_covariance = 10*kf1.observation_covariance,
                      em_vars=['transition_covariance', 'initial_state_covariance'])
    
    kf3 = kf3.em(measurements[:-n_real_time, :], n_iter=5)
    (filtered_state_means, filtered_state_covariances) = kf3.filter(measurements[:-n_real_time,:])
    
    print("Time to build and train kf3: %s seconds" % (time.time() - time_before))
    
    x_now = filtered_state_means[-1, :]
    P_now = filtered_state_covariances[-1, :]
    x_new = np.zeros((n_real_time, filtered_state_means.shape[1]))
    i = 0
    
    for measurement in measurements[-n_real_time:, :]:
        time_before = time.time()
        (x_now, P_now) = kf3.filter_update(filtered_state_mean = x_now,
                                           filtered_state_covariance = P_now,
                                           observation = measurement)
        print("Time to update kf3: %s seconds" % (time.time() - time_before))
        x_new[i, :] = x_now
        i = i + 1
    
    plt.figure(3)
    old_times = range(measurements.shape[0] - n_real_time)
    new_times = range(measurements.shape[0]-n_real_time, measurements.shape[0])
    plt.plot(times, measurements[:, 0], 'bo',
             times, measurements[:, 1], 'ro',
             old_times, filtered_state_means[:, 0], 'b--',
             old_times, filtered_state_means[:, 2], 'r--',
             new_times, x_new[:, 0], 'b-',
             new_times, x_new[:, 2], 'r-')
    
    plt.show()
    

    아래의 플롯은 filter_update 메소드를 사용하여 찾은 3 점을 포함하여 필터 메소드의 성능을 보여줍니다. 도트는 측정 값이고 파선은 필터 학습 기간의 상태 추정치이며 실선은 "온라인"기간의 상태 추정치입니다.

    그리고 타이밍 정보 (내 노트북에서).

    Time to build and train kf3: 0.0677888393402 seconds
    Time to update kf3: 0.00038480758667 seconds
    Time to update kf3: 0.000465154647827 seconds
    Time to update kf3: 0.000463008880615 seconds
    
  2. ==============================

    2.내가 칼만 필터링을 사용하여 볼 수있는 것은 아마도 귀하의 경우에는 올바른 도구가 아닙니다.

    내가 칼만 필터링을 사용하여 볼 수있는 것은 아마도 귀하의 경우에는 올바른 도구가 아닙니다.

    이 방법은 어떨까요? :

    lstInputData = [
        [346, 226 ],
        [346, 211 ],
        [347, 196 ],
        [347, 180 ],
        [350, 2165],  ## noise
        [355, 154 ],
        [359, 138 ],
        [368, 120 ],
        [374, -830],  ## noise
        [346, 90  ],
        [349, 75  ],
        [1420, 67 ],  ## noise
        [357, 64  ],
        [358, 62  ]
    ]
    
    import pandas as pd
    import numpy as np
    df = pd.DataFrame(lstInputData)
    print( df )
    from scipy import stats
    print ( df[(np.abs(stats.zscore(df)) < 1).all(axis=1)] )
    

    여기 출력 :

          0     1
    0    346   226
    1    346   211
    2    347   196
    3    347   180
    4    350  2165
    5    355   154
    6    359   138
    7    368   120
    8    374  -830
    9    346    90
    10   349    75
    11  1420    67
    12   357    64
    13   358    62
          0    1
    0   346  226
    1   346  211
    2   347  196
    3   347  180
    5   355  154
    6   359  138
    7   368  120
    9   346   90
    10  349   75
    12  357   64
    13  358   62
    

    좀 더 자세한 내용은 위의 코드를 참조하십시오.

  3. from https://stackoverflow.com/questions/43377626/how-to-use-kalman-filter-in-python-for-location-data by cc-by-sa and MIT license