Ансамбль фильтра Калмана

Ансамбльный фильтр Калмана ( EnKF ) — это рекурсивный фильтр, подходящий для задач с большим количеством переменных, таких как дискретизация уравнений в частных производных в геофизических моделях. EnKF возник как версия фильтра Калмана для больших задач (по сути, ковариационная матрица заменяется выборочной ковариацией ), и теперь это важный усвоения данных компонент ансамблевого прогнозирования . EnKF связан с фильтром частиц (в этом контексте частица — это то же самое, что и член ансамбля), но EnKF предполагает, что все задействованные распределения вероятностей являются гауссовскими ; когда он применим, он намного более эффективен, чем фильтр твердых частиц .

Введение

Ансамбльный фильтр Калмана (EnKF) представляет собой Монте-Карло реализацию байесовской проблемы обновления : с учетом функции плотности вероятности (PDF) состояния моделируемой системы ( априорная функция , часто называемая прогнозом в науках о Земле) и вероятности данных, Байес Теорема ' используется для получения PDF после того, как была учтена вероятность данных ( апостериорный анализ , часто называемый анализом). Это называется байесовским обновлением. Байесовское обновление сочетается с развитием модели во времени, время от времени включая новые данные. Оригинальный фильтр Калмана , представленный в 1960 году, предполагает, что все PDF-файлы являются гауссовскими (допущение Гаусса), и предоставляет алгебраические формулы для изменения среднего значения и ковариационной матрицы с помощью байесовского обновления, а также формулу для продвижения среднего значения и ковариации во времени при условии, что система линейна. Однако сохранение ковариационной матрицы невозможно с вычислительной точки зрения для многомерных систем. По этой причине были разработаны EnKF. EnKF представляют распределение состояния системы с использованием набора векторов состояния, называемого ансамблем , и заменяют ковариационную матрицу выборочной ковариацией, вычисленной на основе ансамбля. Ансамбль работает так, как будто это случайная выборка , но члены ансамбля на самом деле не являются независимыми , поскольку все они разделяют EnKF. Одним из преимуществ EnKF является то, что продвижение PDF во времени достигается за счет простого продвижения каждого члена ансамбля.

Вывод

Фильтр Калмана

Основная статья: фильтр Калмана

Позволять $\mathbf {x}$ обозначают $n$-мерный вектор состояния модели и предположим, что она имеет гауссово распределение вероятностей со средним значением $\mathbf {\mu }$ и ковариация $Q$, т. е. его PDF

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

Здесь и ниже, $\propto$ означает пропорциональный; PDF-файл всегда масштабируется так, чтобы его интеграл по всему пространству был равен единице. Этот $p(\mathbf {x} )$, называемый предшествующим , со временем изменился в результате запуска модели и теперь должен быть обновлен для учета новых данных. Естественно предположить, что распределение ошибок данных известно; данные должны поступать с оценкой погрешности, иначе они бессмысленны. Здесь данные $\mathbf {d}$ предполагается, что PDF имеет гауссовую форму с ковариацией $R$ и имею в виду $H\mathbf {x}$, где $H$ так называемая матрица наблюдения . Ковариационная матрица $R$ описывает оценку ошибки данных; если случайные ошибки в записях вектора данных $\mathbf {d}$ независимы, $R$ является диагональю, а ее диагональные элементы представляют собой квадраты стандартного отклонения («размера ошибки») ошибки соответствующих элементов вектора данных. $\mathbf {d}$. Значение $H\mathbf {x}$ какова будет ценность данных для государства $\mathbf {x}$ при отсутствии ошибок в данных. Тогда плотность вероятности $p(\mathbf {d} |\mathbf {x} )$ данных $\mathbf {d}$ условие состояния системы $\mathbf {x}$, называемая правдоподобием данных , равна

$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).$

PDF состояния и вероятность данных объединяются, чтобы дать новую плотность вероятности состояния системы. $\mathbf {x}$ зависит от ценности данных $\mathbf {d}$ ( апостериорный ) по теореме Байеса ,

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

Данные $\mathbf {d}$ фиксируется после получения, поэтому обозначим апостериорное состояние через $\mathbf {\hat {x}}$ вместо $\mathbf {x} |\mathbf {d}$ и задний PDF-файл от $p\left(\mathbf {\hat {x}} \right)$. Это можно показать алгебраическими манипуляциями. что апостериорная PDF также является гауссовой,

$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),$

с апостериорным средним значением $\mathbf {\hat {\mu }}$ и ковариация ${\hat {Q}}$ заданный формулами обновления Калмана

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

где

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

– это так называемая матрица усиления Калмана .

Ансамбль Фильтр Калмана

EnKF представляет собой аппроксимацию фильтра Калмана методом Монте-Карло, которая позволяет избежать развития ковариационной матрицы PDF вектора состояния. $\mathbf {x}$. Вместо этого PDF представлен ансамблем

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

$X$ это $n\times N$ матрица, столбцы которой являются членами ансамбля, и называется предшествующим ансамблем . В идеале участники ансамбля должны формировать выборку из предыдущего распределения. Однако члены ансамбля в целом не являются независимыми, за исключением первоначального ансамбля, поскольку каждый шаг EnKF связывает их вместе. Они считаются приблизительно независимыми, и все расчеты ведутся так, как если бы они действительно были независимыми.

Репликация данных $\mathbf {d}$ в $m\times N$ матрица

$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),$

так, чтобы каждый столбец $\mathbf {d} _{i}$ состоит из вектора данных $\mathbf {d}$ плюс случайный вектор из $m$-мерное нормальное распределение $N(0,R)$. Если, кроме того, столбцы $X$ являются выборкой из априорного распределения вероятностей , тогда столбцы

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

сформировать выборку из апостериорного распределения вероятностей . Чтобы увидеть это в скалярном случае с $H=1$: Позволять $x_{i}=\mu +\xi _{i},\;\xi _{i}\sim N(0,\sigma _{x}^{2})$, и $d_{i}=d+\epsilon _{i},\;\epsilon _{i}\sim N(0,\sigma _{d}^{2}).$ Затем

${\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)$.

Первая сумма является апостериорным средним значением, а вторая сумма ввиду независимости имеет дисперсию

$\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}}}$,

что является апостериорной дисперсией.

EnKF теперь получается простой заменой ковариации состояния $Q$ в матрице усиления Калмана $K$ по выборочной ковариации $C$ вычисляется на основе членов ансамбля (так называемая ковариация ансамбля ), то есть: $K=CH^{\mathrm {T} }\left(HCH^{\mathrm {T} }+R\right)^{-1}$

Выполнение

Базовая формулировка

Здесь мы следуем. Предположим, что матрица ансамбля $X$ и матрица данных $D$ такие же, как указано выше. Среднее значение по ансамблю и ковариация равны

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

где

$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},$

и $\mathbf {e}$ обозначает матрицу всех единиц указанного размера.

Задний ансамбль $X^{p}$ затем дается

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

где возмущенная матрица данных $D$ как указано выше.

Обратите внимание, что поскольку $R$ является ковариационной матрицей, она всегда положительно полуопределена и обычно положительно определена , поэтому существует обратное, указанное выше, и формула может быть реализована с помощью разложения Холецкого . В, $R$ заменяется выборочной ковариацией ${\tilde {D}}{\tilde {D}}^{T}/\left(N-1\right)$ где ${\tilde {D}}=D-{\frac {1}{N}}d\,\mathbf {e} _{1\times N}$а обратное заменяется псевдообратным , вычисляемым с использованием разложения по сингулярным значениям (SVD).

Поскольку эти формулы представляют собой матричные операции с доминирующими уровня 3 , операциями они подходят для эффективной реализации с использованием таких программных пакетов, как LAPACK (на компьютерах с последовательной памятью и с общей памятью ) и ScaLAPACK (на компьютерах с распределенной памятью ). Вместо вычисления обратной матрицы и умножения на нее гораздо лучше (в несколько раз дешевле и точнее) вычислить разложение Холецкого матрицы и рассматривать умножение на обратную матрицу как решение линейной системы со многими одновременными правые стороны.

Реализация без матрицы наблюдения

Поскольку мы заменили ковариационную матрицу ансамблевой ковариацией, это приводит к более простой формуле, в которой наблюдения ансамбля используются напрямую без явного указания матрицы $H$. Точнее, определите функцию $h(\mathbf {x} )$ формы

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

Функция $h$ называется функцией наблюдения или, в контексте обратных задач , прямым оператором . Стоимость $h(\mathbf {x} )$ какова будет ценность данных для государства $\mathbf {x}$ при условии, что измерение точное. Тогда задний ансамбль можно переписать в виде

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

где

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

и

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

с

$\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).$

Следовательно, обновление ансамбля можно вычислить путем оценки функции наблюдения $h$ на каждого члена ансамбля один раз и матрица $H$ не обязательно знать явно. Эта формула справедлива и для функции наблюдения $h(\mathbf {x} )=H\mathbf {x+f}$ с фиксированным смещением $\mathbf {f}$, который также не обязательно должен быть известен явно. Приведенная выше формула обычно использовалась для нелинейной функции наблюдения. $h$, например, положение ураганного вихря . В этом случае функция наблюдения по существу аппроксимируется линейной функцией от ее значений в членах ансамбля.

Реализация для большого количества точек данных

Для большого количества $m$ точек данных, умножение на $P^{-1}$ становится узким местом. Следующая альтернативная формула предпочтительна, когда количество точек данных $m$ велика (например, при ассимиляции координатных или пиксельных данных), а ковариационная матрица ошибок данных $R$ является диагональным (это тот случай, когда ошибки данных некоррелированы) или дешевым для разложения (например, полосовым из-за ограниченного ковариационного расстояния). Используя формулу Шермана-Моррисона-Вудбери.

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

с

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

дает

${\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}}$

требующее только решения систем с матрицей $R$ (предполагается, что это дешево) и системы размера $N$ с $m$ правые стороны. Видеть для операций рассчитывается.

Дальнейшие расширения

Описанная здесь версия EnKF предполагает рандомизацию данных. О фильтрах без рандомизации данных см.

Поскольку ковариация ансамбля имеет недостаточный ранг (имеется гораздо больше переменных состояния, обычно миллионы, чем членов ансамбля, обычно менее сотни), она имеет большие члены для пар точек, которые пространственно удалены. Поскольку на самом деле значения физических полей в удаленных местах не так сильно коррелируют , ковариационная матрица искусственно сужается в зависимости от расстояния, что приводит к появлению локализованных алгоритмов EnKF. Эти методы изменяют ковариационную матрицу, используемую в вычислениях, и, следовательно, апостериорный ансамбль больше не состоит только из линейных комбинаций предыдущего ансамбля.

Для нелинейных задач EnKF может создавать апостериорный ансамбль с нефизическими состояниями. Эту проблему можно облегчить с помощью регуляризации , например, штрафования состояний с большими пространственными градиентами .

Для задач с когерентными характеристиками , таких как ураганы , грозы , линии пожара , линии шквалов и фронты дождя , существует необходимость корректировки состояния численной модели путем деформации состояния в пространстве (ее сетки), а также путем аддитивной коррекции амплитуд состояния. . В 2007 году Равела и др. представить совместную модель корректировки положения и амплитуды с использованием ансамблей и систематически выводить последовательное приближение, которое можно применять как к EnKF, так и к другим формулировкам. Их метод не предполагает, что амплитуды и ошибки положения независимы или совместно гауссовы, как это делают другие. Морфинг EnKF использует промежуточные состояния, полученные с помощью методов, заимствованных из регистрации изображений и морфинга , вместо линейных комбинаций состояний.

Формально EnKF полагаются на предположение Гаусса. На практике их также можно использовать для нелинейных задач, где предположение Гаусса может не выполняться. Связанные фильтры, пытающиеся ослабить предположение о Гауссе в EnKF, сохраняя при этом его преимущества, включают фильтры, которые соответствуют PDF состояния с несколькими ядрами Гаусса, фильтры, аппроксимирующие состояние PDF гауссовскими смесями , вариант фильтра частиц с вычислением веса частиц по оценке плотности , и вариант фильтра твердых частиц с с толстыми хвостами PDF-файлом для уменьшения вырождения фильтра частиц .

См. также

• Усвоение данных
• Численный прогноз погоды § Ансамбли
• Фильтр твердых частиц
• Рекурсивная байесовская оценка

Внешние ссылки

• Веб-страница ЭнКФ
• ТОПАЗ, прогнозирование состояния Северного Атлантического океана и арктического морского льда в режиме реального времени с помощью EnKF.
• EnKF-C, компактная платформа для ассимиляции данных в крупномасштабные слоистые геофизические модели с помощью EnKF.
• PDAF - Parallel Data Assimilation Framework - программное обеспечение с открытым исходным кодом для ассимиляции данных, предоставляющее различные варианты EnKF.