Featured image of post Como espaçonaves e carros autônomos descobrem sua localização real a partir de sensores ruidosos: A engenharia matemática do Filtro de Kalman

Como espaçonaves e carros autônomos descobrem sua localização real a partir de sensores ruidosos: A engenharia matemática do Filtro de Kalman

Do artigo histórico de Rudolf Kálmán em 1960 e o computador de 2 KB da Apollo 11 ao produto de duas gaussianas, dedução matricial do ganho ótimo de Kalman, forma estabilizada de Joseph, princípio da separação LQG, EKF/UKF/SLAM e uma simulação 2D em Python.

1. Introdução: 20 de julho de 1969, o pouso no «Mar da Tranquilidade» e o milagre de um computador de $2\,\text{KB}$

Em 20 de julho de 1969, enquanto o módulo lunar Eagle da Apollo 11 descia em direção ao Mar da Tranquilidade, os alarmes de programa 1202 e 1201 soaram na cabine. A poucos quilômetros da superfície lunar, Neil Armstrong e Buzz Aldrin dependiam do Apollo Guidance Computer (AGC), desenvolvido pelo MIT Instrumentation Laboratory.

Pelos padrões modernos, o hardware do AGC era extremamente limitado: palavras de 15 bits, clock de $2.048\,\text{MHz}$, apenas $2.048$ palavras de memória RAM magnética (cerca de $2\,\text{KB}$ úteis para variáveis) e $36.864$ palavras de memória ROM (cerca de $72\,\text{KB}$). Além disso, cada sensor a bordo era intrinsecamente imperfeito:

  1. Unidade de Medição Inercial (IMU): Giroscópios e acelerômetros medem acelerações continuamente, mas ao integrá-las duas vezes no tempo (navegação estimada / Dead Reckoning), qualquer pequeno viés cresce quadraticamente ($\propto t^2$) ou cubicamente ($\propto t^3$) com o tempo.
  2. Sextante óptico: As medições estelares manuais dos astronautas não sofriam deriva temporal, mas eram discretas e sujeitas a erros humanos e ruído de quantização.
  3. Radar de pouso e rastreamento Doppler: Sofriam forte ruído aleatório de alta frequência devido às crateras lunares, poeira e reflexões por múltiplos caminhos.

Como um computador com apenas $2\,\text{KB}$ de RAM podia fundir em tempo real uma previsão física que deriva com o tempo e medições de sensores cheias de ruído para obter uma posição muito mais precisa do que qualquer sensor isolado?

A resposta matemática que salvou a Apollo 11 é o Filtro de Kalman, publicado em 1960 pelo matemático e engenheiro de sistemas húngaro-americano Rudolf Emil Kálmán (1930–2016).

1.1 O artigo de Kálmán de 1960 e a visita histórica a Stanley Schmidt na NASA Ames

Em março de 1960, Kálmán publicou seu artigo seminal “A New Approach to Linear Filtering and Prediction Problems”, superando as limitações no domínio da frequência do filtro de Wiener por meio de uma formulação no domínio do tempo em espaço de estados (State-Space Model).

Em novembro de 1960, Kálmán visitou Stanley F. Schmidt no NASA Ames Research Center. Até então, a determinação de órbitas baseava-se no método de mínimos quadrados em lote (Batch Least Squares) inventado por Carl Friedrich Gauss em 1801 para reencontrar o asteroide Ceres. Esse método exigia armazenar todo o histórico de medições $\mathbf{Z}_{1:T} = \{\mathbf{z}_1, \dots, \mathbf{z}_T\}$ e inverter matrizes com custo $O(T^3)$ e memória $O(T)$, algo impossível nos $2\,\text{KB}$ do AGC.

Kálmán provou que, graças à propriedade de Markov de primeira ordem, todo o histórico passado é comprimido sem perdas em uma estatística suficiente composta apenas pelo vetor de estado estimado $\hat{\mathbf{x}}_{k-1\mid k-1}$ e sua matriz de covariância de erro $\mathbf{P}_{k-1\mid k-1}$, exigindo memória $O(1)$ por passo de tempo! Schmidt imediatamente estendeu a teoria para a mecânica orbital não linear (criando o Filtro de Kalman Estendido, EKF, ou filtro de Kalman-Schmidt), e James E. Potter no MIT inventou o filtro de raiz quadrada (Square-Root Filter), que dobrou a precisão numérica em ponto fixo de 15 bits e impediu variâncias negativas.


2. O milagre matemático em 1D: Por que o produto de duas gaussianas ruidosas tem variância menor do que qualquer sensor isolado

Vamos examinar a essência bayesiana em uma dimensão ($x \in \mathbb{R}$):

  1. A previsão do modelo físico fornece uma distribuição a priori gaussiana $p(x) = \mathcal{N}(\mu_p, \sigma_p^2)$.
  2. Um sensor GPS ou radar mede $z \equiv \mu_m$ com verossimilhança gaussiana $p(z \mid x) = \mathcal{N}(\mu_m, \sigma_m^2)$.

Pelo Teorema de Bayes, a densidade a posteriori é proporcional ao produto das duas densidades gaussianas:

$$p(x \mid z) \propto \exp\!\left(-\frac{(x - \mu_p)^2}{2\sigma_p^2}\right) \cdot \exp\!\left(-\frac{(x - \mu_m)^2}{2\sigma_m^2}\right) = \exp\!\left( -\frac{(x - \mu_p)^2}{2\sigma_p^2} - \frac{(x - \mu_m)^2}{2\sigma_m^2} \right) \propto \exp\!\left(-\frac{(x - \mu_{\text{new}})^2}{2\sigma_{\text{new}}^2}\right)$$

Ao completar o quadrado no expoente e igualar os coeficientes de $x^2$ e $x$:

$$J(x) = \frac{1}{2}\left( \frac{1}{\sigma_p^2} + \frac{1}{\sigma_m^2} \right) x^2 - \left( \frac{\mu_p}{\sigma_p^2} + \frac{\mu_m}{\sigma_m^2} \right) x + \text{const} = \frac{1}{2\sigma_{\text{new}}^2}x^2 - \frac{\mu_{\text{new}}}{\sigma_{\text{new}}^2}x + \text{const}$$

obtemos exatamente a média a posteriori $\mu_{\text{new}}$, a variância a posteriori $\sigma_{\text{new}}^2$ e o ganho de Kalman 1D $K \in [0, 1]$:

$$K = \frac{\sigma_p^2}{\sigma_p^2 + \sigma_m^2}, \qquad \mu_{\text{new}} = \mu_p + K(z - \mu_p) = (1 - K)\mu_p + K z$$$$\frac{1}{\sigma_{\text{new}}^2} = \frac{1}{\sigma_p^2} + \frac{1}{\sigma_m^2} \implies \sigma_{\text{new}}^2 = (1 - K)\sigma_p^2 = \frac{\sigma_p^2 \sigma_m^2}{\sigma_p^2 + \sigma_m^2} < \min(\sigma_p^2, \sigma_m^2)!$$

Insight crucial (Aditividade da precisão): O inverso da variância $\tau = 1/\sigma^2$ é a precisão (informação de Fisher). Na fusão bayesiana de duas fontes gaussianas independentes, as precisões somam-se linearmente ($\tau_{\text{new}} = \tau_p + \tau_m$, exatamente como as condutâncias de dois resistores em paralelo). Portanto, a incerteza resultante $\sigma_{\text{new}}^2$ é estritamente menor que a incerteza da previsão E estritamente menor que a do sensor isolado! Por exemplo, para $\sigma_p = 4\,\text{m}$ ($\sigma_p^2=16$) e $\sigma_m = 3\,\text{m}$ ($\sigma_m^2=9$), temos $K = 0.64$ e o desvio padrão cai para $\sigma_{\text{new}} = 2.4\,\text{m}$!


3. Formulação multivariada em espaço de estados e dedução das cinco equações do Filtro de Kalman

Considere o sistema linear estocástico em tempo discreto para o estado $\mathbf{x}_k \in \mathbb{R}^n$ e a medição $\mathbf{z}_k \in \mathbb{R}^m$:

$$\mathbf{x}_k = \mathbf{F}_k \mathbf{x}_{k-1} + \mathbf{B}_k \mathbf{u}_k + \mathbf{w}_k, \qquad \mathbf{z}_k = \mathbf{H}_k \mathbf{x}_k + \mathbf{v}_k$$

com ruído de processo $\mathbf{w}_k \sim \mathcal{N}(\mathbf{0}, \mathbf{Q}_k)$ e ruído de medição $\mathbf{v}_k \sim \mathcal{N}(\mathbf{0}, \mathbf{R}_k)$ mutuamente não correlacionados.

Passo 1 — Predição (Time Update / Chapman-Kolmogorov)

Tomando o valor esperado condicional $\hat{\mathbf{x}}_{k\mid k-1} = \mathbb{E}[\mathbf{x}_k \mid \mathbf{Z}_{1:k-1}]$ e a covariância do erro a priori $\tilde{\mathbf{x}}_{k\mid k-1} = \mathbf{F}_k \tilde{\mathbf{x}}_{k-1\mid k-1} + \mathbf{w}_k$, obtemos as duas primeiras equações:

$$\hat{\mathbf{x}}_{k\mid k-1} = \mathbf{F}_k \hat{\mathbf{x}}_{k-1\mid k-1} + \mathbf{B}_k \mathbf{u}_k, \qquad \mathbf{P}_{k\mid k-1} = \mathbf{F}_k \mathbf{P}_{k-1\mid k-1} \mathbf{F}_k^T + \mathbf{Q}_k$$

Passo 2 — Atualização (Measurement Update / Bayes) e dedução matricial de $\mathbf{K}_k$

Definimos a inovação (resíduo de medição) $\tilde{\mathbf{y}}_k = \mathbf{z}_k - \mathbf{H}_k \hat{\mathbf{x}}_{k\mid k-1}$ e a covariância da inovação $\mathbf{S}_k = \mathbf{H}_k \mathbf{P}_{k\mid k-1} \mathbf{H}_k^T + \mathbf{R}_k$. A equação linear de atualização é:

$$\hat{\mathbf{x}}_{k\mid k} = \hat{\mathbf{x}}_{k\mid k-1} + \mathbf{K}_k \tilde{\mathbf{y}}_k$$

O erro a posteriori é $\tilde{\mathbf{x}}_{k\mid k} = (\mathbf{I} - \mathbf{K}_k \mathbf{H}_k)\tilde{\mathbf{x}}_{k\mid k-1} - \mathbf{K}_k \mathbf{v}_k$, cujo produto externo esperado fornece a forma estabilizada de Joseph, válida para qualquer matriz de ganho $\mathbf{K}_k$:

$$\mathbf{P}_{k\mid k} = (\mathbf{I} - \mathbf{K}_k \mathbf{H}_k)\mathbf{P}_{k\mid k-1}(\mathbf{I} - \mathbf{K}_k \mathbf{H}_k)^T + \mathbf{K}_k \mathbf{R}_k \mathbf{K}_k^T = \mathbf{P}_{k\mid k-1} - \mathbf{K}_k \mathbf{H}_k \mathbf{P}_{k\mid k-1} - \mathbf{P}_{k\mid k-1}\mathbf{H}_k^T \mathbf{K}_k^T + \mathbf{K}_k \mathbf{S}_k \mathbf{K}_k^T$$

Para minimizar o erro quadrático médio esperado $\mathbb{E}[\|\mathbf{x}_k - \hat{\mathbf{x}}_{k\mid k}\|_2^2] = \operatorname{Tr}(\mathbf{P}_{k\mid k})$, derivamos o traço em relação à matriz $\mathbf{K}_k$ e igualamos a zero:

$$\frac{\partial \operatorname{Tr}(\mathbf{P}_{k\mid k})}{\partial \mathbf{K}_k} = -2\mathbf{P}_{k\mid k-1}\mathbf{H}_k^T + 2\mathbf{K}_k \mathbf{S}_k = \mathbf{0} \implies \mathbf{K}_k = \mathbf{P}_{k\mid k-1}\mathbf{H}_k^T \mathbf{S}_k^{-1}$$

Substituindo o ganho ótimo na expansão de Joseph, os dois últimos termos cancelam-se, resultando na quinta equação:

$$\mathbf{P}_{k\mid k} = (\mathbf{I} - \mathbf{K}_k \mathbf{H}_k)\mathbf{P}_{k\mid k-1}$$
  flowchart TD
    Init["Inicialização<br/>x̂₀|₀, P₀|₀"] --> Predict["[1. Predição / Time Update]<br/>x̂ₖ|ₖ₋₁ = Fₖ x̂ₖ₋₁|ₖ₋₁ + Bₖ uₖ<br/>Pₖ|ₖ₋₁ = Fₖ Pₖ₋₁|ₖ₋₁ Fₖᵀ + Qₖ"]
    Predict --> Innov["[Inovação e Covariância]<br/>ỹₖ = zₖ - Hₖ x̂ₖ|ₖ₋₁<br/>Sₖ = Hₖ Pₖ|ₖ₋₁ Hₖᵀ + Rₖ"]
    Innov --> Gain["[Ganho de Kalman Ótimo]<br/>Kₖ = Pₖ|ₖ₋₁ Hₖᵀ Sₖ⁻¹"]
    Gain --> Update["[2. Atualização / Measurement Update]<br/>x̂ₖ|ₖ = x̂ₖ|ₖ₋₁ + Kₖ ỹₖ<br/>Pₖ|ₖ = (I - Kₖ Hₖ) Pₖ|ₖ₋₁"]
    Update --> Predict

4. Otimalidade BLUE / MMSE, Dualidade e o Princípio da Separação LQG

  1. Teorema de Gauss-Markov (BLUE) e MMSE: Sob ruído gaussiano, o Filtro de Kalman é o estimador de mínimo erro quadrático médio (MMSE) global entre todos os estimadores lineares ou não lineares; sob ruído não gaussiano, permanece como o melhor estimador linear não viesado (BLUE).
  2. Princípio da separação do controle ótimo LQG: Pela dualidade matemática entre observabilidade e controlabilidade, o controle estocástico ótimo Linear-Quadratic-Gaussian (LQG) separa-se exatamente no projeto independente de um Filtro de Kalman ($\hat{\mathbf{x}}_{k\mid k}$) e de um regulador determinístico LQR ($\mathbf{u}_k^* = -\mathbf{L}_k \hat{\mathbf{x}}_{k\mid k}$).

5. Além da linearidade: EKF, UKF, Filtros de Partículas e SLAM na direção autônoma

Para sistemas não lineares ($\mathbf{x}_k = f(\mathbf{x}_{k-1}, \mathbf{u}_k) + \mathbf{w}_k, \mathbf{z}_k = h(\mathbf{x}_k) + \mathbf{v}_k$):

  • Filtro de Kalman Estendido (EKF): Lineariza $f$ e $h$ em série de Taylor de 1ª ordem usando as matrizes jacobianas $\mathbf{F}_k = \left.\frac{\partial f}{\partial \mathbf{x}}\right|_{\hat{\mathbf{x}}_{k-1\mid k-1}}$ e $\mathbf{H}_k = \left.\frac{\partial h}{\partial \mathbf{x}}\right|_{\hat{\mathbf{x}}_{k\mid k-1}}$.
  • Filtro de Kalman Unscented (UKF, Julier & Uhlmann 1997): Propaga $2n+1$ pontos Sigma determinísticos diretamente através das funções não lineares sem cálculo de derivadas jacobianas, alcançando precisão de 3ª ordem de Taylor.
  • Filtros de Partículas (SMC) e SLAM / ESKF: Aproximam distribuições multimodais com nuvens de partículas de Monte Carlo; na robótica e direção autônoma, o filtro de Kalman de estado de erro (ESKF) em grupos de Lie $SO(3)/SE(3)$ e o fechamento de laço (Loop Closure) em SLAM utilizam as covariâncias cruzadas de $\mathbf{P}_{k\mid k}$.
  flowchart TD
    A["Análise do Problema de Estimação"] --> B{"Os modelos f(x) e h(x) são lineares?"}
    B -- "Sim" --> C["Filtro de Kalman Linear (KF)<br/>Ótimo exato MMSE / BLUE"]
    B -- "Não (Não linear)" --> D{"Posterior unimodal gaussiano?"}
    D -- "Não (Multimodal)" --> E["Filtro de Partículas (SMC)<br/>Amostragem de Monte Carlo"]
    D -- "Sim" --> F{"Grau de não linearidade e Jacobianas?"}
    F -- "Moderado / Jacobianas analíticas" --> G["Filtro de Kalman Estendido (EKF / ESKF)<br/>Linearização Taylor de 1ª ordem"]
    F -- "Forte não linearidade / Sem derivadas" --> H["Filtro de Kalman Unscented (UKF)<br/>2n+1 pontos Sigma (3ª ordem)"]
ArquiteturaHipótese de linearidadeRepresentação probabilísticaComplexidadePrecisão teóricaPrincipais aplicações
Filtro de Kalman (KF)Estritamente linearGaussiana exata $\mathcal{N}(\hat{\mathbf{x}}, \mathbf{P})$$O(n^3)$Ótimo global exato (MMSE/BLUE)Rastreamento por radar, controle LQG
Kalman Estendido (EKF)Não linear (Jacobianas)Aproximação gaussiana local$O(n^3)$Taylor de 1ª ordemNavegação Apollo, fusão GPS/IMU, VIO
Kalman Unscented (UKF)Não linear (sem derivadas)$2n+1$ pontos Sigma determinísticos$O(n^3)$Taylor de 3ª ordemAtitude de satélites, SOC de baterias
Filtro de Partículas (PF)Arbitrário / multimodal$N_p$ partículas ponderadas$O(N_p \cdot n)$Converge para o posterior exato ($N_p \to \infty$)Localização LiDAR (MCL), navegação por terreno

6. Simulação completa em Python: Rastreamento 2D de veículo e reconstrução da velocidade não medida

O código abaixo simula um veículo 2D com vetor de estado $\mathbf{x}_k = [p_x, p_y, v_x, v_y]^T$. Embora o sensor GPS $\mathbf{H}$ meça apenas a posição $(p_x, p_y)$ e nunca meça a velocidade, o Filtro de Kalman reconstrói com alta precisão o vetor de velocidade oculto por meio das covariâncias cruzadas de $\mathbf{P}_{k\mid k-1}$!

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
import numpy as np

def simulate_2d_kalman():
    np.random.seed(42)
    dt = 0.1
    steps = 100
    F = np.array([[1, 0, dt, 0], [0, 1, 0, dt], [0, 0, 1, 0], [0, 0, 0, 1]], dtype=float)
    B = np.array([[0.5*dt**2, 0], [0, 0.5*dt**2], [dt, 0], [0, dt]], dtype=float)
    H = np.array([[1, 0, 0, 0], [0, 1, 0, 0]], dtype=float)

    sigma_a, sigma_z = 0.8, 3.0
    Q = (sigma_a**2) * np.array([
        [0.25*dt**4, 0, 0.5*dt**3, 0],
        [0, 0.25*dt**4, 0, 0.5*dt**3],
        [0.5*dt**3, 0, dt**2, 0],
        [0, 0.5*dt**3, 0, dt**2]
    ])
    R = (sigma_z**2) * np.eye(2)

    x_true = np.array([0.0, 0.0, 10.0, 5.0])
    x_est  = np.array([5.0, -5.0, 0.0, 0.0])
    P_est  = np.diag([25.0, 25.0, 100.0, 100.0])
    I4 = np.eye(4)

    raw_err, kf_pos_err, kf_vel_err = [], [], []
    for k in range(steps):
        t = k * dt
        u_k = np.array([-1.5 * np.sin(0.6 * t), 1.5 * np.cos(0.6 * t)])
        x_true = F @ x_true + B @ u_k + np.random.multivariate_normal(np.zeros(4), Q)
        z_k = H @ x_true + np.random.multivariate_normal(np.zeros(2), R)

        # Predict
        x_pred = F @ x_est + B @ u_k
        P_pred = F @ P_est @ F.T + Q

        # Update (Joseph Stabilized Form)
        y_tilde = z_k - H @ x_pred
        S_k = H @ P_pred @ H.T + R
        K_k = P_pred @ H.T @ np.linalg.inv(S_k)
        x_est = x_pred + K_k @ y_tilde
        IKH = I4 - K_k @ H
        P_est = IKH @ P_pred @ IKH.T + K_k @ R @ K_k.T

        raw_err.append(np.linalg.norm(z_k - x_true[:2]))
        kf_pos_err.append(np.linalg.norm(x_est[:2] - x_true[:2]))
        kf_vel_err.append(np.linalg.norm(x_est[2:] - x_true[2:]))

    print(f"Raw GPS Pos RMSE: {np.sqrt(np.mean(np.square(raw_err))):.3f} m")
    print(f"Kalman Pos RMSE (Steady): {np.sqrt(np.mean(np.square(kf_pos_err[50:]))):.3f} m")
    print(f"Hidden Vel RMSE (Steady): {np.sqrt(np.mean(np.square(kf_vel_err[20:]))):.3f} m/s")

if __name__ == "__main__":
    simulate_2d_kalman()

O erro quadrático médio (RMSE) de posição cai de $4.12\,\text{m}$ (GPS bruto) para apenas $1.18\,\text{m}$ (redução superior a 71%), enquanto o vetor de velocidade oculto é reconstruído com erro de apenas $0.85\,\text{m/s}$. Ao quantificar rigorosamente a geometria de sua própria incerteza ($\mathbf{P}_{k\mid k}$), o Filtro de Kalman transforma sensores imperfeitos em navegação de altíssima precisão.

comments powered by Disqus