Filtro de Kalman de conjunto

Dejar${\displaystyle \mathbf {x} }$denotan el${\displaystyle n}$vector de estado de dimensión de un modelo, y supongamos que tiene una distribución de probabilidad gaussiana con media${\displaystyle \mathbf {\mu } }$y covarianza${\displaystyle Q}$, es decir, su PDF es

${\displaystyle p(\mathbf {x} )\propto \exp \left(-{\frac {1}{2}}(\mathbf {x} -\mathbf {\mu } )^{\mathrm {T} }Q^{-1}(\mathbf {x} -\mathbf {\mu } )\right).}$

Aquí y abajo,${\displaystyle \propto }$significa proporcional; una PDF siempre se escala de manera que su integral sobre todo el espacio sea uno.${\displaystyle p(\mathbf {x} )}$, llamado prior , evolucionó con el tiempo al ejecutar el modelo y ahora debe actualizarse para tener en cuenta los nuevos datos. Es natural suponer que se conoce la distribución de errores de los datos; los datos deben venir con una estimación de error, de lo contrario no tienen sentido. Aquí, los datos${\displaystyle \mathbf {d} }$Se supone que tiene una función de densidad de probabilidad gaussiana con covarianza${\displaystyle R}$y significa${\displaystyle H\mathbf {x} }$, dónde${\displaystyle H}$es la llamada matriz de observación . La matriz de covarianza${\displaystyle R}$describe la estimación del error de los datos; si los errores aleatorios en las entradas del vector de datos${\displaystyle \mathbf {d} }$son independientes,${\displaystyle R}$es diagonal y sus entradas diagonales son los cuadrados de la desviación estándar ("tamaño del error") del error de las entradas correspondientes del vector de datos.${\displaystyle \mathbf {d} }$. El valor${\displaystyle H\mathbf {x} }$es cuál sería el valor de los datos para el estado${\displaystyle \mathbf {x} }$en ausencia de errores en los datos. Entonces la densidad de probabilidad${\displaystyle p(\mathbf {d} |\mathbf {x} )}$de los datos${\displaystyle \mathbf {d} }$condicional del estado del sistema${\displaystyle \mathbf {x} }$, llamado probabilidad de datos , es

${\displaystyle p\left(\mathbf {d} |\mathbf {x} \right)\propto \exp \left(-{\frac {1}{2}}(\mathbf {d} -H\mathbf {x} )^{\mathrm {T} }R^{-1}(\mathbf {d} -H\mathbf {x} )\right).}$

La función de densidad de probabilidad del estado y la verosimilitud de los datos se combinan para obtener la nueva densidad de probabilidad del estado del sistema.${\displaystyle \mathbf {x} }$condicionado al valor de los datos${\displaystyle \mathbf {d} }$(la posterior ) por el teorema de Bayes ,

${\displaystyle p\left(\mathbf {x} |\mathbf {d} \right)\propto p\left(\mathbf {d} |\mathbf {x} \right)p(\mathbf {x} ).}$

Los datos${\displaystyle \mathbf {d} }$queda fijo una vez recibido, por lo que denotamos el estado posterior por${\displaystyle \mathbf {\hat {x}} }$en lugar de${\displaystyle \mathbf {x} |\mathbf {d} }$y la PDF posterior por${\displaystyle p\left(\mathbf {\hat {x}} \right)}$. Se puede demostrar mediante manipulaciones algebraicas ^{[ 5 ]} que la PDF posterior también es gaussiana,

${\displaystyle p\left(\mathbf {\hat {x}} \right)\propto \exp \left(-{\frac {1}{2}}(\mathbf {\hat {x}} -\mathbf {\hat {\mu }} )^{\mathrm {T} }{\hat {Q}}^{-1}(\mathbf {\hat {x}} -\mathbf {\hat {\mu }} )\right),}$

con la media posterior${\displaystyle \mathbf {\hat {\mu }} }$y covarianza${\displaystyle {\hat {Q}}}$dadas por las fórmulas de actualización de Kalman

${\displaystyle \mathbf {\hat {\mu }} =\mathbf {\mu } +K\left(\mathbf {d} -H\mathbf {\mu } \right),\quad {\hat {Q}}=\left(I-KH\right)Q,}$

dónde

${\displaystyle K=QH^{\mathrm {T} }\left(HQH^{\mathrm {T} }+R\right)^{-1}}$

es la denominada matriz de ganancia de Kalman .

El EnKF es una aproximación de Monte Carlo del filtro de Kalman, que evita la evolución de la matriz de covarianza de la función de densidad de probabilidad del vector de estado.${\displaystyle \mathbf {x} }$En cambio, la PDF está representada por un conjunto

${\displaystyle X=\left[\mathbf {x} _{1},\ldots ,\mathbf {x} _{N}\right]=\left[\mathbf {x} _{i}\right].}$

${\displaystyle X}$es un${\displaystyle n\times N}$La matriz cuyas columnas son los miembros del conjunto se denomina conjunto previo . Idealmente, los miembros del conjunto formarían una muestra de la distribución previa. Sin embargo, los miembros del conjunto no son independientes en general , excepto en el conjunto inicial, ya que cada paso del EnKF los vincula. Se consideran aproximadamente independientes, y todos los cálculos se realizan como si realmente lo fueran.

Replicar los datos${\displaystyle \mathbf {d} }$en un${\displaystyle m\times N}$matriz

${\displaystyle D=\left[\mathbf {d} _{1},\ldots ,\mathbf {d} _{N}\right]=\left[\mathbf {d} _{i}\right],\quad \mathbf {d} _{i}=\mathbf {d} +\mathbf {\epsilon _{i}} ,\quad \mathbf {\epsilon _{i}} \sim N(0,R),}$

de modo que cada columna${\displaystyle \mathbf {d} _{i}}$consta del vector de datos${\displaystyle \mathbf {d} }$más un vector aleatorio del${\displaystyle m}$distribución normal dimensional${\displaystyle N(0,R)}$. Si, además, las columnas de${\displaystyle X}$son una muestra de la distribución de probabilidad previa , luego las columnas de

${\displaystyle {\hat {X}}=X+K(D-HX)}$

formar una muestra de la distribución de probabilidad posterior . Para ver esto en el caso escalar con${\displaystyle H=1}$: Dejar${\displaystyle x_{i}=\mu +\xi _{i},\;\xi _{i}\sim N(0,\sigma _{x}^{2})}$, y${\displaystyle d_{i}=d+\epsilon _{i},\;\epsilon _{i}\sim N(0,\sigma _{d}^{2}).}$ Entonces

${\displaystyle {\hat {x}}_{i}=\left({\frac {1/\sigma _{x}^{2}}{1/\sigma _{x}^{2}+1/\sigma _{d}^{2}}}\mu +{\frac {1/\sigma _{d}^{2}}{1/\sigma _{x}^{2}+1/\sigma _{d}^{2}}}d\right)+\left({\frac {1/\sigma _{x}^{2}}{1/\sigma _{x}^{2}+1/\sigma _{d}^{2}}}\xi _{i}+{\frac {1/\sigma _{d}^{2}}{1/\sigma _{x}^{2}+1/\sigma _{d}^{2}}}\epsilon _{i}\right)}$.

La primera suma es la media posterior, y la segunda suma, en vista de la independencia, tiene una varianza.

${\displaystyle \left({\frac {1/\sigma _{x}^{2}}{1/\sigma _{x}^{2}+1/\sigma _{d}^{2}}}\right)^{2}\sigma _{x}^{2}+\left({\frac {1/\sigma _{d}^{2}}{1/\sigma _{x}^{2}+1/\sigma _{d}^{2}}}\right)^{2}\sigma _{d}^{2}={\frac {1}{1/\sigma _{x}^{2}+1/\sigma _{d}^{2}}}}$,

que es la varianza posterior.

El EnKF ahora se obtiene simplemente reemplazando la covarianza del estado.${\displaystyle Q}$en la matriz de ganancia de Kalman${\displaystyle K}$por la covarianza de la muestra${\displaystyle C}$calculada a partir de los miembros del conjunto (llamada covarianza del conjunto ), ^{[ 6 ]} es decir:${\displaystyle K=CH^{\mathrm {T} }\left(HCH^{\mathrm {T} }+R\right)^{-1}}$

Aquí seguimos. ^{[ 7 ] [ 8 ]} Supongamos que la matriz de conjunto${\displaystyle X}$y la matriz de datos${\displaystyle D}$son como se indicó anteriormente. La media del conjunto y la covarianza son

${\displaystyle E\left(X\right)={\frac {1}{N}}\sum _{k=1}^{N}\mathbf {x} _{k},\quad C={\frac {AA^{T}}{N-1}},}$

dónde

${\displaystyle A=X-E\left(X\right)\mathbf {e} _{1\times N}=X-{\frac {1}{N}}\left(X\mathbf {e} _{N\times 1}\right)\mathbf {e} _{1\times N},}$

y${\displaystyle \mathbf {e} }$denota la matriz de todos los unos del tamaño indicado.

El conjunto posterior${\displaystyle X^{p}}$entonces se da por

${\displaystyle X^{p}=X+CH^{T}\left(HCH^{T}+R\right)^{-1}(D-HX),}$

donde la matriz de datos perturbada${\displaystyle D}$es como se indicó anteriormente.

Tenga en cuenta que desde${\displaystyle R}$es una matriz de covarianza, siempre es semidefinida positiva y usualmente definida positiva , por lo que la inversa anterior existe y la fórmula puede implementarse mediante la descomposición de Cholesky . ^{[ 9 ]} En, ^{[ 7 ] [ 8 ]}${\displaystyle R}$se reemplaza por la covarianza de la muestra${\displaystyle {\tilde {D}}{\tilde {D}}^{T}/\left(N-1\right)}$dónde${\displaystyle {\tilde {D}}=D-{\frac {1}{N}}d\,\mathbf {e} _{1\times N}}$y la inversa se reemplaza por una pseudoinversa , calculada utilizando la descomposición en valores singulares (SVD).

Dado que estas fórmulas son operaciones matriciales con operaciones de Nivel 3 dominantes , ^{[ 10 ]} son ​​adecuadas para una implementación eficiente utilizando paquetes de software como LAPACK (en computadoras de memoria serial y compartida ) y ScaLAPACK (en computadoras de memoria distribuida ). ^{[ 9 ]} En lugar de calcular la inversa de una matriz y multiplicarla por ella, es mucho mejor (varias veces más barato y también más preciso) calcular la descomposición de Cholesky de la matriz y tratar la multiplicación por la inversa como la solución de un sistema lineal con muchos segundos miembros simultáneos. ^{[ 10 ]}

Dado que hemos reemplazado la matriz de covarianza por la covarianza de conjunto, esto conduce a una fórmula más simple donde las observaciones de conjunto se utilizan directamente sin especificar explícitamente la matriz.${\displaystyle H}$. Más específicamente, defina una función${\displaystyle h(\mathbf {x} )}$de la forma

${\displaystyle h(\mathbf {x} )=H\mathbf {x} .}$

La función${\displaystyle h}$se denomina función de observación o, en el contexto de los problemas inversos , operador directo . El valor de${\displaystyle h(\mathbf {x} )}$es cuál sería el valor de los datos para el estado${\displaystyle \mathbf {x} }$suponiendo que la medición sea exacta. Entonces el conjunto posterior se puede reescribir como

${\displaystyle X^{p}=X+{\frac {1}{N-1}}A\left(HA\right)^{T}P^{-1}(D-HX)}$

dónde

${\displaystyle HA=HX-{\frac {1}{N}}\left(\left(HX\right)\mathbf {e} _{N\times 1}\right)\mathbf {e} _{1\times N},}$

y

${\displaystyle P={\frac {1}{N-1}}HA\left(HA\right)^{T}+R,}$

con

${\displaystyle \left[HA\right]_{i}=H\mathbf {x} _{i}-H{\frac {1}{N}}\sum _{j=1}^{N}\mathbf {x} _{j}\ =h\left(\mathbf {x} _{i}\right)-{\frac {1}{N}}\sum _{j=1}^{N}h\left(\mathbf {x} _{j}\right).}$

En consecuencia, la actualización del conjunto se puede calcular evaluando la función de observación.${\displaystyle h}$en cada miembro del conjunto una vez y la matriz${\displaystyle H}$No es necesario conocerlo explícitamente. Esta fórmula también es válida ^{[ 9 ]} para una función de observación.${\displaystyle h(\mathbf {x} )=H\mathbf {x+f} }$con un desplazamiento fijo${\displaystyle \mathbf {f} }$, que tampoco necesita conocerse explícitamente. La fórmula anterior se ha utilizado comúnmente para una función de observación no lineal.${\displaystyle h}$, como la posición de un vórtice de huracán . ^{[ 11 ]} En ese caso, la función de observación se aproxima esencialmente mediante una función lineal a partir de sus valores en los miembros del conjunto.

Para un gran número${\displaystyle m}$de puntos de datos, la multiplicación por${\displaystyle P^{-1}}$se convierte en un cuello de botella. La siguiente fórmula alternativa es ventajosa cuando el número de puntos de datos${\displaystyle m}$es grande (como cuando se asimilan datos en cuadrícula o de píxeles) y la matriz de covarianza de errores de los datos${\displaystyle R}$es diagonal (que es el caso cuando los errores de los datos no están correlacionados), o barato de descomponer (como en bandas debido a la distancia de covarianza limitada). Usando la fórmula de Sherman-Morrison-Woodbury ^{[ 12 ]}

${\displaystyle (R+UV^{T})^{-1}=R^{-1}-R^{-1}U(I+V^{T}R^{-1}U)^{-1}V^{T}R^{-1},}$

con

${\displaystyle U={\frac {1}{N-1}}HA,\quad V=HA,}$

da

${\displaystyle {\begin{aligned}P^{-1}&=\left(R+{\frac {1}{N-1}}HA\left(HA\right)^{T}\right)^{-1}\ =\\&=R^{-1}\left[I-{\frac {1}{N-1}}\left(HA\right)\left(I+\left(HA\right)^{T}R^{-1}{\frac {1}{N-1}}\left(HA\right)\right)^{-1}\left(HA\right)^{T}R^{-1}\right],\end{aligned}}}$

lo cual requiere únicamente la solución de sistemas con la matriz${\displaystyle R}$(se supone que es barato) y de un sistema de tamaño${\displaystyle N}$con${\displaystyle m}$Lados derechos. Véase ^{[ 9 ]} para el recuento de operaciones.