Известия РАН. Теория и системы управления, 2020, № 1, стр. 54-67

Введение. Основная проблема оптимальной нелинейной фильтрации марковских сигналов состоит в большой сложности реализации в необходимом темпе со временем абсолютно оптимального фильтра (АОФ) [1–5], устанавливающего явную зависимость оценки от всех накопленных измерений. Он является бесконечномерной динамической (рекуррентной) системой с распределенными параметрами, так как его состояние определяется функцией апостериорной плотности вероятности. Поэтому осуществляющий непосредственную потраекторную реализацию АОФ современный фильтр частиц [6–8] требует применения весьма мощного вычислителя из-за использования им громоздкого последовательного метода Монте-Карло. На практике это заставляет применять субоптимальные конечномерные алгоритмы оценивания вроде обобщенного фильтра Калмана, который представляет собой простейшее линеаризованное приближение к АОФ. Однако точность таких приближений может оказаться недостаточной, в частности на большом отрезке времени из-за рекуррентного накопления методических и вычислительных погрешностей [9, 10]. Необходимость же вычислять кроме вектора самой оценки еще и матрицу его ковариаций делает даже такие, лишь ковариационные приближения к АОФ, достаточно сложными для реализации, особенно в случае значительного числа оцениваемых переменных [11].

Альтернативой медленному АОФ являются конечномерные фильтры оптимальной структуры (ФОС) различных порядков, представляющие собой обобщения условно оптимального фильтра заданной структуры и фиксированного порядка [3, 12, 13]. Непрерывно-дискретные версии ФОС для марковских непрерывных (диффузионных) сигналов приведены в [14, 15], а непрерывный ФОС для кусочно-непрерывных (диффузионно-скачкообразных) сигналов описан в [16]. При конечномерном подходе задача нелинейной фильтрации разбивается на две части, и наиболее сложная из них, состоящая в получении уравнений фильтра, решается заранее, до появления измерений, на этапе проектирования фильтра. В результате реализация ФОС, состоящая в рекуррентном получении им оценок по поступающим измерениям, требует лишь решения (моделирования) этих уравнений. Трудоемкость этой операции определяется числом уравнений и видом их нелинейностей.

Однако и перечисленные непрерывно-дискретные ФОС обладают рядом недостатков. Точность имеющего непрерывный прогноз фильтра малого порядка (ФМП) [14] ограничена из-за неизменной размерности его вектора состояния, равной числу оцениваемых переменных. Порядок фильтра с конечной памятью (ФКП) и с кусочно-постоянным прогнозом [15] уже можно выбирать, так как он кратен размерности вектора измерений, но оценка этого фильтра зависит лишь от нескольких последних измерений.

В настоящей работе предлагается процедура построения достаточно простого для реализации рекуррентного фильтра большого порядка (ФБП) с кусочно-полиномиальным прогнозом. Его время памяти, в отличие от ФКП, уже бесконечно, а противоречие между точностью и сложностью можно регулировать выбором порядка фильтра. Вектор оценки по-прежнему ищется как наилучшая функция от последнего измерения и вектора состояния фильтра. Но последний теперь формируется из векторов нескольких предыдущих оценок, увеличивая этим, по сравнению с ФМП, оперативную память фильтра, что позволяет получать лучшую точность за счет расширения допустимого множества оценок. При этом старые измерения не забываются, ибо информация о них аккумулируется в оценках. Кроме того, класс диффузионных оцениваемых сигналов расширяется на диффузионно-скачкообразные.

Показано, что синтез нового ФБП сводится к нахождению заранее соответствующей условной плотности вероятности из рекуррентной цепочки преобразований типа прогноз–коррекция. Ее основу составляют интегродифференциальное уравнение Колмогорова–Феллера в частных производных (прогноз) и интегральная формула Байеса (коррекция). Иначе синтезировать ФБП можно численно методом Монте-Карло, но с громоздким построением на каждом измерительном такте гистограммы искомой структурной вектор-функции фильтра, зависящей от большого числа аргументов. Поэтому рассмотрено и построение двух традиционных численно-аналитических ковариационных приближений к ФБП. Демонстрируется пример оценивания состояния стохастического осциллятора Ван-дер-Поля.

1. Постановка задачи непрерывно-дискретной фильтрации. Рассмотрим следующую скрытую марковскую модель системы наблюдения.

Пусть изменяющийся во времени $t \in [0,\,T]$ кусочно-непрерывный n-мерный оцениваемый полезный сигнал ${{X}_{t}} \in {{\mathbb{R}}^{n}}$, непрерывный справа ${{X}_{t}} = {{X}_{{t + }}}$ и имеющей предел слева ${{X}_{{t - }}}$, является вектором состояния своего формирующего фильтра. Последний, представляющий собой объект наблюдения, возмущается не только непрерывными гауссовскими, но и импульсными пуассоновскими белыми шумами. Он описывается стохастическим дифференциальным уравнением Ито, которое в интегральной форме записи имеет следующий вид:

Здесь X₀ – случайное начальное значение сигнала, характеризуемое плотностью вероятности ${{p}_{0}}(x)$; $a(t,x)$ – вектор-функция сноса; $B(t,x)$ – матричная функция диффузии; W_t – стандартный винеровский процесс; ${{S}_{t}}({{X}_{{t - }}}) = {{X}_{t}} - {{X}_{{t - }}}$ – случайный вектор амплитуды скачка случайного процесса X_t в момент времени t, который определяется условной плотностью вероятности $\xi (t{\text{,}}s{\kern 1pt} {\text{|}}{\kern 1pt} x)$; ${{\Theta }_{1}},{{\Theta }_{2}}, \ldots $ – пуассоновский поток моментов времени скачков, которые зависят от текущего значения сигнала и имеют условную интенсивность $\mu (t,x)$:

$P_{0}^{t}(x)$ – считающий число этих скачков за промежуток времени [0, t) целочисленный условный процесс с законом распределения Пуассона

Отметим, что в уравнении (1.1) интеграл по времени понимается как среднеквадратический, а по винеровскому процессу ${{W}_{t}}$ – как стохастический интеграл Ито. Известно, что плотность вероятности $p(t,x)$ кусочно-непрерывного сигнала ${{X}_{t}}$ удовлетворяет интегродифференциальному уравнению Колмогорова–Феллера [2, 17, 18], которое является обобщением более известного дифференциального уравнения Фоккера–Планка–Колмогорова на класс диффузионно-скачкообразных процессов.

Пусть также в известные тактовые моменты времени

производятся мгновенные безынерционные измерения вектора ${{X}_{t}}$, в том числе неполные или неточные, определяемые по формуле измерителя

Здесь ${{Y}_{k}} \in {{\mathbb{R}}^{m}}$ – m-мерный вектор измерения, ${{c}_{k}}(x,v)$ – вектор-функция измерителя, V_k – вектор независимого дискретного белого шума с плотностью вероятности ${{q}_{k}}(v)$.

Относительно соотношений (1.1), (1.2) сделаем следующие естественные предположения.

Предположение 1. Уравнение (1.1) имеет сильное решение ${{X}_{t}}$ и оно единственно. Достаточные для этого условия приведены в [17, 18].

Предположение 2. Случайные процесс ${{X}_{t}}$ и его измерения ${{Y}_{k}}$ имеют конечные вторые моменты

где ${\rm M}$ – оператор математического ожидания, ${{\left\| \varepsilon \right\|}^{2}} = {{\varepsilon }^{{\text{T}}}}\varepsilon $ – квадрат евклидовой нормы.

Предположение 3. Для простоты изложения будем считать все случайные величины абсолютно-непрерывными, что позволяет характеризовать их плотностями вероятности.

Требуется в каждый момент времени из любого полуинтервала между двумя соседними тактовыми моментами времени найти оценку ${{\hat {X}}_{t}}$ вектора ${{X}_{t}}$ как функцию всех накопленных к этому времени измерений:

оптимальную в смысле минимума среднеквадратической ошибки оценивания

Здесь ${{C}_{t}} = C_{t}^{{\text{T}}} \succ 0$ – симметрическая положительно определенная матрица весовых коэффициентов. Отметим, что критерий (1.4) обеспечивает и несмещенность оценки

что гарантирует отсутствие в ней систематической погрешности. Поэтому этот критерий служит для минимизации разброса оценки относительно оцениваемой величины.

2. Выбор уравнения нового фильтра. Скорректируем приведенную постановку задачи, для чего сначала опишем известный результат ее решения.

2.1. Построение абсолютно оптимального фильтра. Если на класс оценок (1.3), определяемых функциями ${{\hat {x}}_{{t,k}}}( \cdot )$ растущего со временем числа аргументов $m(k + 1)$, никаких ограничений не накладывать, то оптимальной в смысле (1.4) является функция апостериорного среднего:

Здесь ${{\rho }_{k}}( \cdot )$ – апостериорная плотность вероятности, а интегралы здесь и далее берутся по всему евклидову пространству соответствующей размерности:

При этом сама плотность ${{\rho }_{k}}( \cdot )$ от предыдущих оценивающих функций (2.1) и соответствующих им оценок ${{\hat {X}}_{t}}$ не зависит и на каждом промежутке $({{t}_{k}},{{t}_{{k + 1}}}]$ между измерениями удовлетворяет интегродифференциальному уравнению Колмогорова–Феллера в частных производных (уравнение прогноза). В момент очередного измерения ${{t}_{{k + 1}}}$ конечное сечение ${{\rho }_{k}}(t_{{k + 1}}^{ - },x{\kern 1pt} {\text{|}}{\kern 1pt} y_{0}^{k})$ решения этого уравнения по формуле Байеса (уравнение коррекции) пересчитывается в начальное для следующего промежутка $({{t}_{{k + 1}}},{{t}_{{k + 2}}}]$ сечение ${{\rho }_{{k + 1}}}(t_{{k + 1}}^{ + },x{\kern 1pt} {\text{|}}{\kern 1pt} y_{0}^{{k + 1}})$. Реализовать указанные вычисления в темпе с поступлением измерений, как правило, затруднительно.

Поэтому откажемся от требования явной зависимости оценки от всех измерений (1.3), которая и приводит к указанной выше сложности АОФ. Отметим, что поиск подобной (1.3) явной зависимости тактовых оценок ${{\hat {X}}_{{{{t}_{k}}}}}$ диффузионного сигнала X_t лишь от нескольких последних измерений ${{Y}_{k}},{{Y}_{{k - 1}}}, \ldots ,{{Y}_{{k - l}}}$ привел к существенно более простому, хотя и менее точному, непрерывно-дискретному ФОС [15], однако время его памяти конечно.

2.2. Фильтр с полиномиальным прогнозом. С целью получить столь же простой дискретный ФОС, но уже с бесконечным временем памяти, предлагается на каждом измерительном такте k искать наилучшую оценку ${{Z}_{k}}$ только тактового значения ${{X}_{{t_{k}^{{\, - }}}}}$ n-мерного диффузионно-скачкообразного вектора ${{X}_{t}}$. Для возможной экономии вычислительных затрат можно также оценивать лишь некоторую его наиболее интересную $n{\kern 1pt} '$-мерную часть $X_{t}^{'} \in {{\mathbb{R}}^{{n'}}}$, которую составим, без ограничения общности, из первых $n{\kern 1pt} ' \leqslant n$ компонент вектора ${{X}_{t}}$. Поэтому критерий (1.4) оптимальности оценки всюду, в любой момент времени, заменим подобным, но более слабым, тактовым условием оптимальности:

Саму же оценку ${{Z}_{k}}$ будем искать в виде ее явной зависимости от последнего измерения ${{Y}_{k}}$ и от не более чем $l \in \mathbb{N}$ предыдущих оценок $Z_{{k - l}}^{{k - 1}} = \left[ {{{Z}_{{k - 1}}},{{Z}_{{k - 2}}}, \ldots ,{{Z}_{{k - l}}}} \right]$, а именно

определяя наилучшую функцию ${{f}_{k}}( \cdot )$ из условия (2.2).

Прогноз $\hat {X}_{t}^{'}$ информационного вектора $X_{t}^{'}$ между измерениями построим с помощью одной или нескольких оптимальных тактовых оценок ${{Z}_{k}}$ эвристически, например в виде экстраполирующего полинома. Так, прогноз можно сделать постоянным, линейным и т.д.:

Таким образом, для оценивания случайного процесса (1.1) по дискретным измерениям (1.2) предлагается синтезировать непрерывно-дискретный фильтр с оптимальной тактовой оценкой (2.3) и полиномиальным прогнозом (2.4). В отличие от ФКП [15] здесь старые измерения не забываются, ибо информация о них аккумулируется в оценках, поэтому время памяти нового фильтра как динамической системы бесконечно. Подобный фильтр, но для дискретных марковских сигналов, описан в [19].

2.3. Фильтр с кусочно-полиномиальным прогнозом. Для получения более точного прогноза на каждом отрезке времени $[{{t}_{k}},{{t}_{{k + 1}}}]$ между соседними измерениями введем аналогично [15] дополнительно L межтактовых точек $\tau _{k}^{i}$_, $i = \overline {1,L} $, обозначив границы отрезка как $\tau _{k}^{0}$ и $\tau _{k}^{{L + 1}}$:

В этих точках найдем и межтактовые оценки $Z_{k}^{i}$ как их явные зависимости от l тактовых оценок:

определяя функции $g_{k}^{i}( \cdot )$ из аналогичного (2.2) условия оптимальности:

Используя эти дополнительные оценки $Z_{k}^{i}$, можно строить какой-либо из полиномиальных прогнозов типа (2.4) и на каждом межтактовом промежутке отдельно:

3. Вектор состояния и порядок предлагаемого фильтра. Уравнение (2.3), связывающее текущую оценку Z_k с предыдущими, является рекуррентной формулой (разностным уравнением) высокого порядка. Преобразуем его в систему уравнений первого порядка, для чего соберем l последних оценок в блочный вектор-столбец растущей вначале размерности:

Процесс изменения этого вектора ${{U}_{k}}$ учитываемых оценок можно записать рекуррентно:

с помощью следующей вектор-функции накопления или обновления оценок:

Здесь С – матрица удаления из блочного вектора-столбца ${{U}_{{k - 1}}}$ его последнего, устаревшего, блока ${{Z}_{{k - l}}}$, тогда как Е, О – единичная и нулевая матрицы соответственно. Тогда основное уравнение фильтра (2.3) и дополнительные выражения для межтактовых оценок (2.5) принимают вид

Замечание 1. Вместо (3.4) можно искать и подобные (2.3) непосредственные зависимости межтактовых оценок $Z_{k}^{i}$ и от последнего измерения ${{Y}_{k}}$, заменив (2.5) соотношениями

С учетом (3.1) эти формулы будут аналогичны (3.3):

В итоге исходные уравнения фильтра типа вход-выход (2.3), (2.5) представлены эквивалентными соотношениями (3.2)–(3.4) типа вход-состояние-выход. Действительно, разностное уравнение первого порядка (3.2) представляет собой уравнение состояния фильтра, а соотношения (3.3), (3.4) – формулы его выходов. При этом само состояние фильтра определяется блочным вектором (3.1), а его максимальная размерность $p = ln{\kern 1pt} '$, не  изменяющаяся  с l-го такта $k = l - 1$, является порядком фильтра. За счет увеличения коэффициента кратности l этот порядок может быть сделан сколь угодно большим. Оптимизации же по критерию (2.2) подлежат только функции выхода фильтра ${{f}_{k}}( \cdot )$, $f_{k}^{i}( \cdot )$, тогда как функция его состояния ${{s}_{k}}( \cdot )$ фиксирована. В результате выше доказано следующее утверждение.

Лемма. Определяемый рекуррентными формулами (2.3), (2.5) фильтр является конечномерным с фиксированным уравнением состояния (3.2) и с оптимизируемыми формулами выходов (3.3), (3.4), а его вектор состояния (3.1) состоит из желаемого числа l запоминаемых им оценок, так что порядок фильтра равен $ln{\kern 1pt} '$, где $n{\kern 1pt} '$ – размерность вектора оценки.

Замечание 2. Из-за неизменности функции состояния ${{s}_{k}}( \cdot )$ предлагаемый алгоритм оценивания относится к типу полуоптимальных ФОС произвольного порядка [20]. От ФКП новый фильтр отличается принципиальной заменой прошлых измерений $Y_{{\max (0,k - l)}}^{{k - 1}}$ на прошлые оценки $Z_{{\max (0,k - l)}}^{{k - 1}}$. В частном же случае, при коэффициенте кратности l = 1, когда уравнения (2.3), (2.5) принимают наиболее простой вид

этот фильтр вырождается в известный ФМП [15, 21].

Поэтому дадим следующее определение.

Определение. Рекуррентную зависимость тактовых оценок от измерений (2.3), как и эквивалентную ей пару из уравнения состояния (3.2) и формулы выхода (3.3), назовем фильтром большого порядка $p = ln{\kern 1pt} '$ с полиномиальным прогнозом (2.4), а эти же соотношения, но с дополнительными формулами для межтактовых прогнозов (2.5) или (3.4) будем именовать фильтром с кусочно-полиномиальным прогнозом (2.7).

4. Нахождение функций выхода фильтра. Подставляя искомые формулы для оценок (2.3), (2.5) в соответствующие им критерии оптимальности (2.2), (2.6), из известной теоремы о наилучшей среднеквадратической регрессии получаем следующий результат, подобный (2.1).

Теорема 1 (об оптимальной структуре ФБП). При выполнении предположений 1–3 наилучшие среди всех измеримых по Борелю функций выхода фильтра существуют и находятся как следующие условные средние:

а соответствующие им оценки ${{Z}_{k}}$, $Z_{k}^{i}$ являются несмещенными.

Однако эти соотношения лишь выражают искомые функции через три условных плотности вероятности: начальную ${{\rho }_{0}}({{x}_{0}}{\kern 1pt} {\text{|}}{\kern 1pt} {{y}_{0}})$, корректирующую ${{\rho }_{k}}(t_{k}^{ - },x{\kern 1pt} {\text{|}}{\kern 1pt} {{y}_{k}},{{u}_{{k - 1}}})$ и прогнозирующую ${{\pi }_{k}}(t,x{\kern 1pt} {\text{|}}{\kern 1pt} {{u}_{k}})$, $t \in [{{t}_{k}},{{t}_{{k + 1}}})$. Используя формулу Байеса, первую из них найдем по известной плотности ${{p}_{0}}(x)$ начального состояния объекта (1.1):

а каждую корректирующую плотность ${{\rho }_{k}}({{x}_{{{{t}_{k}}}}}{\kern 1pt} {\text{|}}{\kern 1pt} {{y}_{k}},{{u}_{{k - 1}}})$ представим через конечное, при $t = t_{k}^{ - }$, сечение предыдущей прогнозирующей ${{\pi }_{{k - 1}}}(t,x{\kern 1pt} {\text{|}}{\kern 1pt} {{u}_{{k - 1}}})$, которая справедлива лишь на промежутке $t \in [{{t}_{{k - 1}}},{{t}_{k}})$:

В двух последних соотношениях ${{\beta }_{k}}({{y}_{k}}|{{x}_{k}})$ – получаемая из формулы измерителя (1.2) функция правдоподобия:

где $\delta ( \cdot )$ – функция Дирака, а символом $numerator$ здесь и далее обозначен числитель соответствующей дроби. Аналогично и прогнозирующую плотность ${{\pi }_{k}}(t,x{\kern 1pt} {\text{|}}{\kern 1pt} {{u}_{k}})$ выразим через совместную плотности вероятности ${{r}_{k}}(t,x,{{u}_{k}})$ случайных состояний объекта наблюдения X_t и фильтра U_k:

Наконец, используя марковское свойство определяемого уравнением Ито (1.1) диффузионно-скачкообразного процесса X_t, можно показать, подобно [15], что каждая из плотностей вероятности ${{r}_{k}}( \cdot )$ на соответствующем ей промежутке времени между измерениями удовлетворяет известному уравнению Колмогорова–Феллера

Здесь K_x, F_x – два прямых производящих оператора процесса: дифференциальный оператор Фоккера–Планка–Колмогорова

с вектором ${{\nabla }_{x}}$ градиента по х и оператором следа матрицы tr, а также интегральный оператор Феллера

Отметим, что при недифференцируемости по переменной х функции сноса $a(t,x)$ один раз, а функции диффузии $B(t,x)$ дважды, решение уравнения (4.8) понимают в обобщенном смысле Бубнова–Галеркина [22, 23].

Из приведенных в (3.2), (3.3) связей между случайными переменными ${{U}_{0}} = {{Z}_{0}}$, ${{Z}_{0}} = {{f}_{0}}({{Y}_{0}})$ и известных свойств плотностей вероятности следует, что начальным условием для уравнения (4.8) на первом промежутке $t \in [{{t}_{0}},{{t}_{1}})$ является плотность вероятности

Здесь функция начальной оценки ${{f}_{0}}({{y}_{0}})$ уже известна из (4.2), (4.4).

Аналогично находится и начальное условие ${{r}_{k}}({{t}_{k}},x,{{u}_{k}})$ для уравнения (4.8) на каждом из последующих межтактовых промежутков времени $t \in [{{t}_{k}},{{t}_{{k + 1}}})$, $k \geqslant 1$. Оно оказывается связанным с конечным значением ${{r}_{{k - 1}}}(t_{k}^{ - },x,{{u}_{{k - 1}}})$ решения этого уравнения на предыдущем промежутке $t \in [{{t}_{{k - 1}}},{{t}_{k}})$ из-за изменения в тактовый момент времени $t = {{t}_{k}}$ вектора состояния фильтра с ${{U}_{{k - 1}}}$ на ${{U}_{k}}$ по формуле (3.2). Действительно, в результате появления нового измерения Y_k и получения по (3.3) очередной оценки Z_k с помощью уже найденной функции ${{f}_{k}}({{y}_{k}},{{u}_{{k - 1}}})$ последняя добавляется к вектору ${{U}_{{k - 1}}}$ с возможным удалением из него на этапе обновления уже устаревшей оценки ${{Z}_{{k - l}}}$. Так как на этапе накопления оценок ${{U}_{k}} = Z_{0}^{k}$, а на этапе их обновления ${{U}_{k}} = Z_{{k - l + 1}}^{k}$, то и сечения плотностей вероятности ${{r}_{k}}( \cdot )$, ${{r}_{{k - 1}}}( \cdot )$ в один и тот же момент времени t_k оказываются связанными одним из двух следующих соотношений:

Здесь функция ${{f}_{k}}( \cdot )$ уже известна из (4.2), (4.5), (4.7), а интегрирование по переменной ${{z}_{{k - l}}}$ во втором из этих выражений соответствует удалению из вектора ${{U}_{{k - 1}}} = Z_{{k - l}}^{{k - 1}}$ его блока ${{Z}_{{k - l}}}$.

5. Алгоритм точного синтеза фильтра. Приведенные выше соотношения позволяют определить следующую последовательность действий. Они выполняются заранее, до появления результатов измерений, а потому могут быть осуществлены на достаточно мощном вычислителе. Последнее выгодно отличает этот фильтр, как и другие ФОС, от классического АОФ.

Первым способом нахождения оптимальных тактовых функций фильтра (4.2) является чередование процесса решения уравнения Колмогорова–Феллера (4.8) на очередном межтактовом промежутке времени $t \in [{{t}_{k}},{{t}_{{k + 1}}})$, $k = 0,1, \ldots $, с подходящим тактовым пересчетом (4.10) конечного значения полученной плотности вероятности ${{r}_{k}}(t_{{k + 1}}^{ - },x,{{u}_{k}})$ в начальное условие ${{r}_{{k + 1}}}({{t}_{{k + 1}}},x,{{u}_{{k + 1}}})$ для решения этого же уравнения на следующем промежутке. Для выполнения такого пересчета требуется тактовая функция ${{f}_{{k + 1}}}({{y}_{{k + 1}}},{{u}_{k}})$, поэтому каждое полученное на k-м промежутке решение ${{r}_{k}}( \cdot )$ этого уравнения по формуле (4.7) преобразуется в прогнозирующую плотность вероятности ${{\pi }_{k}}( \cdot )$, конечное сечение которой по (4.5) пересчитывается в корректирующую плотность ${{\rho }_{{k + 1}}}( \cdot )$. Последняя и позволяет по (4.2) найти функцию ${{f}_{{k + 1}}}( \cdot )$, которая, во-первых, является результатом синтеза фильтра на данном такте и, во-вторых, позволяет осуществить упомянутый тактовый пересчет ${{r}_{k}}( \cdot )$ в ${{r}_{{k + 1}}}( \cdot )$ по (4.10). Начальным же условием для решения уравнения (4.8) на промежутке $t \in [{{t}_{0}},{{t}_{1}})$ является сразу получаемая по (4.1), (4.4), (4.9) плотность вероятности ${{r}_{0}}( \cdot )$. Одновременно с этим каждое межтактовое сечение прогнозирующей плотности ${{\pi }_{k}}( \cdot )$ в момент времени $\tau _{k}^{i}$ позволяет по (4.3) найти еще и функцию $g_{k}^{i}({{u}_{k}})$ вычисления дополнительной оценки $Z_{k}^{i}$, служащей для построения более точного кусочно-полиномиального прогноза (2.7).

Вторым численным способом синтеза фильтра является получение практически точных, вследствие выполнения этих расчетов заранее, оптимальных функций условных средних (4.1)–(4.3) с помощью применяемого пошагово во времени метода Монте-Карло. Для этого необходимо на каждом измерительном такте k осуществить многократное статистическое моделирование уравнений объекта (1.1), измерителя (1.2) и фильтра (3.2), (3.3), с целью вычисления достаточно больших пакетов реализаций случайных величин $\{ {{X}_{{t_{k}^{{\, - }}}}}\} $, $\{ {{Y}_{k}}\} $, $\{ {{U}_{{k - 1}}}\} $. При этом интегрирование уравнения Ито можно выполнить с помощью одной из известных стохастических разностных схем типа Рунге–Кутты, например методом Эйлера–Маруямы. Затем следует выполнить k-ю обработку результатов моделирования, находя гистограмму оптимальной функции условного среднего ${{f}_{k}}({{y}_{k}},{{u}_{{k - 1}}})$ и осуществляя ее сглаживание каким-либо из известных способов. Знание этой функции позволяет продолжить процесс моделирования фильтра на следующий такт. При необходимости одновременно с этим следует фиксировать и межтактовые значения сигнала $\{ {{X}_{{\tau _{k}^{i} - }}}\} $ для аналогичного получения дополнительных функций $g_{k}^{i}({{u}_{k}})$.

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

6. Гауссовское приближение к фильтру. Используем приведенную в [15] процедуру построения этого известного приближения к ФКП, который, как было отмечено выше, отличается от предлагаемого ФБП только видом вектора состояния U_k.

Аппроксимируем числители дробей формул Байеса (4.4), (4.5), (4.7), а именно одну условную и две безусловные плотность вероятности:

нормальной плотностью $N( \cdot \parallel m,D)$ с соответствующими им условными или безусловными параметрами типа среднего m и ковариации D. Тогда по свойству гауссовского закона распределения получим, что определяемые этими дробями условные плотности ${{\rho }_{k}}(t_{k}^{{\, - }},x{\kern 1pt} {\text{|}}{\kern 1pt} {{y}_{k}},{{u}_{{k - 1}}})$, ${{\pi }_{k}}(t,x{\kern 1pt} {\text{|}}{\kern 1pt} {{u}_{k}})$, ${{\rho }_{0}}({{x}_{0}}{\kern 1pt} {\text{|}}{\kern 1pt} {{y}_{0}})$ тоже почти гауссовские, а их параметры, среди которых присутствуют и искомые условные средние (4.1)–(4.3), выражаются через m, D по теореме о нормальной корреляции. Поэтому оптимальные формулы выхода фильтра (3.3) получают приближения, лишь линейные по измерениям ${{Y}_{k}}$, ${{Y}_{0}}$. Кроме того, благодаря следующей связи между плотностями и виду (4.6) функции правдоподобия ${{\beta }_{k}}({{y}_{k}}{\kern 1pt} {\text{|}}{\kern 1pt} x)$, три из пяти параметров гауссовской аппроксимации условной плотности ${{\eta }_{k}}(t_{k}^{ - },x,{{y}_{k}}{\kern 1pt} {\text{|}}{\kern 1pt} {{u}_{{k - 1}}})$ удается выразить через два других ее параметра с помощью функции ${{c}_{k}}(x,v)$ измерителя (1.2). В результате аналогично [15] получим такой алгоритм вычисления тактовых и межтактовых оценок.

Теорема 2 (об уравнениях гауссовского ФБП). Если у двух условных средних измерителя (1.2)

существуют гауссовские моменты (характеристики их статистической линеаризации по Казакову) то имеют место такие уравнения субоптимального фильтра:

Здесь ${{\Lambda }_{k}}$ – вектор тактового прогноза, $\Lambda _{k}^{'},\,\,T_{k}^{'}$ – матрицы из первых $n{\kern 1pt} '$ строк $n$-строчных матриц ${{\Lambda }_{k}}$, ${{T}_{k}}$, $ \oplus $ – символ псевдообращения матрицы по Муру–Пенроузу, ${{H}_{0}},\,\,{{e}_{0}}$ – начальные параметры фильтра, ${{\Gamma }_{k}},\,\,{{\kappa }_{k}},\,\,{{T}_{k}}$ – тактовые параметры, $\Gamma _{k}^{i},\,\,\kappa _{k}^{i}$ – межтактовые параметры. При этом все эти параметры выражаются через два первых момента m, D известных случайных величин по формулам:

Отметим, что эти параметры тоже определяются методом Монте-Карло путем обработки результатов многократного потактового статистического моделирования уравнений объекта (1.1), измерителя (1.2) и фильтра (6.2), но теперь уже легко, без построения гистограмм.

7. Линеаризованное приближение к фильтру. К сожалению, необходимая для построения гауссовского фильтра интегральная процедура нахождения характеристик статистической линеаризации нелинейностей (6.1) довольно сложна, а сами они в реальных задачах оценивания существуют далеко не всегда [24]. Поэтому применим и другое известное приближение, существенно более простое, а потому и более популярное. Аналогично [15] легко получим такое утверждение.

Следствие (о функциях коррекции линеаризованного ФБП). Если нелинейность ${{c}_{k}}(x,v)$ измерителя (1.2) дифференцируема по обоим аргументам и для нее в точке $({{\Lambda }_{k}},m_{k}^{v})$, где ${{\Lambda }_{k}}$ – прогноз вектора ${{X}_{{t_{k}^{ - }}}}$, $m_{k}^{v}$ – среднее значение помехи V_k, справедлива формула Тейлора

то функции гауссовской коррекции (6.1) аппроксимируются следующими выражениями:

Здесь $C_{k}^{x}(x) = {{\left. {{{\partial {{c}_{k}}(x,v)} \mathord{\left/ {\vphantom {{\partial {{c}_{k}}(x,v)} {\partial x}}} \right. \kern-0em} {\partial x}}} \right|}_{{v = m_{k}^{v}}}}$, $C_{k}^{v}(x) = {{\left. {{{\partial {{c}_{k}}(x,v)} \mathord{\left/ {\vphantom {{\partial {{c}_{k}}(x,v)} {\partial v}}} \right. \kern-0em} {\partial v}}} \right|}_{{v = m_{k}^{v}}}}$ – сечения матриц Якоби частных производных, ${{R}_{k}} = {\text{cov}}\,[{{V}_{k}},{{V}_{k}}]$ – матрица ковариации помехи измерения.

В результате уравнения линеаризованного ФБП также имеют вид (6.2), но с легко получаемыми функциями (7.1), а его параметры по-прежнему вычисляются по формулам (6.3).

Замечание 3. Недостатки двух рассмотренных численно-аналитических приближений к ФОС как линеаризованного (дифференцируемость нелинейностей и малость отклонения от точки линеаризации), так и гауссовского (сложность аналитического получения гауссовских моментов нелинейностей и возможность их несуществования) можно преодолеть применением занимающих промежуточное положение между ними чисто численных приближений типа кубатурного [25] и нечувствительного (unscented) или, иначе, сигма-точечного [9, 26, 27].

8. Сравнение нового и классического ковариационных фильтров. Для выяснения принципиальных отличий субоптимальных ФБП (6.2) и ФКП [15] от аналогичного приближения к непрерывно-дискретному АОФ приведем уравнения последнего. При этом ограничимся для простоты рассмотрением задачи оценивания чисто диффузионного сигнала, полагая, что скачки у процесса X_t отсутствуют, так что ${{X}_{{t - }}} = {{X}_{t}}$. Более того, пусть ФБП оценивает весь вектор X_t, т.е. $n{\kern 1pt} ' = n$, вследствие чего его оценки ${{Z}_{k}}$, $Z_{k}^{i}$ можно обозначить как ${{\hat {X}}_{{{{t}_{k}}}}}$, $\hat {X}_{{{{t}_{k}}}}^{i}$.

Линеаризованная версия АОФ хорошо известна как обобщенный (extended) фильтр Калмана (ОФК), а гауссовская – как менее известный фильтр нормальной аппроксимации (ФНА). Их уравнения также различаются лишь видом структурных функций. Любой из них строит прогноз между измерениями $t \in [{{t}_{k}},{{t}_{{k + 1}}})$ непрерывно, интегрируя автономную (замкнутую) систему обыкновенных дифференциальных уравнений для оценки ${{\hat {X}}_{t}} = {\rm M}[{{X}_{t}}{\kern 1pt} {\text{|}}{\kern 1pt} Y_{0}^{k}]$ и матрицы ${{P}_{t}} = \operatorname{cov} [{{X}_{t}}{\kern 1pt} {\text{|}}{\kern 1pt} Y_{0}^{k}]$ апостериорных ковариаций оцениваемого вектора [3, с. 354]:

Затем в каждый тактовый момент времени ${{t}_{k}}$ конечные значения ${{\hat {X}}_{{t_{k}^{ - }}}}$, ${{P}_{{t_{k}^{ - }}}}$ решений системы (8.1) на предыдущем промежутке $t \in [{{t}_{{k - 1}}},{{t}_{k}})$ корректируются измерением Y_k в начальные условия ${{\hat {X}}_{{{{t}_{k}}}}}$, ${{P}_{{{{t}_{k}}}}}$ для следующего промежутка $t \in [{{t}_{k}},{{t}_{{k + 1}}})$ по формулам [3, с. 458]:

Здесь первоначальные значения являются детерминированными:

Таким образом, кроме тех же, что и в ФБП, функций коррекции ${{h}_{k}}( \cdot ),\,\,{{G}_{k}}( \cdot ),$ ${{F}_{k}}( \cdot )$, теперь необходимы еще и три функции прогноза $\tau \,( \cdot )$, $A\,( \cdot )$, $\Theta \,( \cdot )$. В гауссовском приближении для ФНА они тоже находятся как подобные (6.1) характеристики статистической линеаризации, но уже от функций сноса и диффузии уравнения объекта (1.1):

а в линеаризованном приближении для ОФК они так же, как и (7.1), имеют более простой вид:

Однако систему дифференциальных уравнений (8.1), второе из которых является уравнением типа Риккати, приходится решать численно с достаточно малым шагом, обеспечивая при этом симметричность и положительную определенность матрицы ковариаций, что совсем не просто [9, 10], а n-мерный вектор оценки ${{\hat {X}}_{t}}$ и различные элементы $n \times n{\text{ - }}$матрицы ковариации ${{P}_{{{{t}_{k}}}}}$ образуют довольно большой вектор состояния фильтра, так что его порядок равен $n(n + 3){\text{/}}2$. Причем здесь сэкономить на оценивании лишь части переменных состояния не удается.

Сопоставляя уравнения нового фильтра (6.2) с известными (8.1), (8.2), можно заметить следующие принципиальные преимущества приближенных ФОС.

1. Тактовые оценки субоптимальных фильтров вычисляются с помощью одной и той же (общей для них) корректирующей функции

но по разным значениям ее аргументов:

Здесь тактовый АОФ-прогноз ${{\hat {X}}_{{t_{k}^{ - }}}}$, полученный интегрированием системы уравнений (8.1), заменен вектором ФОС-прогноза Λ_k, который легко находится по вектору состояния этого фильтра ${{U}_{{k - 1}}}~ = \hat {X}_{{{{t}_{{\max (0,k - l)}}}}}^{{{{t}_{{k - 1}}}}}$ как ${{\Lambda }_{k}} = {{\Gamma }_{k}}{{U}_{{k - 1}}} + {{\kappa }_{k}}$ с помощью заранее вычисленных параметров Γ_k, ${{\kappa }_{k}}$, а потому является несмещенным. Кроме того, получаемая из той же системы (8.1) тактовая матрица ${{P}_{{t_{k}^{ - }}}} = {\text{cov}}[{{X}_{{t_{k}^{ - }}}}{\kern 1pt} {\text{|}}{\kern 1pt} Y_{0}^{{k - 1}}]$ заменена в ФОС заранее найденной детерминированной матрицей T_k = = ${\text{cov}}[{{X}_{{t_{k}^{ - }}}} - {{\Lambda }_{k}}]$, благодаря чему сложное уравнение типа Риккати решать тоже не нужно.

2. На промежутке времени между измерениями прогноз АОФ строится решением сложной системы дифференциальных уравнений (8.1), тогда как прогноз ФБП определяется простой экстраполяцией (2.4) или (2.7) оценок тактовых ${{\hat {X}}_{{{{t}_{k}}}}}$ или межтактовых $\hat {X}_{{{{t}_{k}}}}^{i}$. Последние вычисляются как ${\hat {X}_{{{{t}_{k}}}}^{i} = \Gamma _{k}^{i}{{U}_{k}} + \kappa _{k}^{i}}$ и являются несмещенными.

Проведенный анализ позволяет сделать следующий вывод.

Утверждение (о сравнении приближений к ФБП и к АОФ). Приближенным ФОС не требуется интегрировать систему уравнений приближенных АОФ для вектора прогноза и матрицы его ковариаций (8.1), а потому не нужно находить и три функции прогноза по (8.3) или по (8.4). Вместо этого необходимо заранее найти по (6.3) и использовать тактовые параметры ${{\Gamma }_{k}},{{\kappa }_{k}},{{T}_{k}}$ и межтактовые параметры $\Gamma _{k}^{i},\,\,\kappa _{k}^{i}$, которые компенсируют отсутствие информации о ковариации ошибки оценивания P_t.

9. Пример численного сравнения субоптимальных фильтров²². Для анализа точности приближений к АОФ и к ФОС рассмотрим задачу оценивания двумерного вектора состояния стохастической версии нелинейного осциллятора Ван-дер-Поля. Соответствующая система уравнений Ито (1.1), которая здесь приведена в нестрогой форме Ланжевена, имеет вид

где ${{\dot {W}}_{t}}$ – гауссовский белый шум. На рис. 1 представлена фазовая плоскость (${{X}_{{1,t}}}$, ${{X}_{{2,t}}}$) с двумя предельными циклами этого осциллятора, возмущенным (при ${{\dot {W}}_{t}} \ne 0$) и невозмущенным (при ${{\dot {W}}_{t}}$ = 0), которые образуют соответствующие петли гистерезиса.

Пусть обе переменные состояния этого осциллятора измеряются с аддитивными погрешностями, так что дискретный измеритель (1.2) является двумерным и нелинейным по X_1,t:

Точность сравниваемых фильтров определялась методом Монте-Карло с усреднением по выборке в 800 реализаций. Статистическое моделирование дифференциальных уравнений осциллятора и фильтра с непрерывным прогнозом проводилось методом типа Рунге–Кутты на отрезке [0, 3] с шагом интегрирования Δt_инт = 0.01, интервал времени между измерениями принят равным Δt_изм = 0.5, а между кусочно-постоянными дискретными прогнозами Δt_прогн = 0.1. Структурные функции линеаризованных и гауссовских фильтров для этого примера найдены в [15].

На рис. 2 приведены графики среднеквадратических отклонений (СКО) ошибок Se1 оценивания переменной X_1,t несколькими линеаризованными (Л) непрерывно-дискретными фильтрами с разными видами прогнозов. Сравниваются четыре фильтра:

1) Л-АОФ (ОФК) (8.1), (8.2), (8.4), прогноз которого с необходимостью непрерывен,

2) Л-ФМП с непрерывным прогнозом [14],

3) Л-ФБП кратности l = 1 с кусочно-постоянными прогнозами (2.7) при $L = 4$,

4) Л-ФКП кратности l = 4 [15] с такими же кусочно-постоянными прогнозами.

Из рис. 2 видно, что по точности Л-АОФ, который в этом примере имеет порядок p = = $2 \cdot (2 + 3){\text{/}}2$ = 5, со временем все более проигрывает обоим Л-ФМП порядка p = 2 и тем более 4-кратному Л-ФКП порядка $р = 4 \cdot 2 = 8$.

Время же реализации новых фильтров по сравнению с классическим Л-АОФ тоже оказалось лучше:

у Л-ФМП с непрерывным прогнозом оно меньше в 1.6 раза, так как последний, в отличие от Л-АОФ, не интегрирует три дифференциальных уравнения для ковариаций,

у Л-ФМП и Л-ФКП с кусочно-постоянными прогнозами оно меньше даже в 2.8 раза, ибо они не интегрируют никакие дифференциальные уравнения, даже для прогноза, причем большой порядок Л-ФКП на скорость оценивания практически не повлиял.

Таким образом, этот пример не только подтверждает сформулированный в разд. 8 вывод о вычислительных преимуществах линеаризованных ФБП и ФКП, но и демонстрирует даже некоторое их превосходство по точности перед ОФК.

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

Получены формулы для оптимальных функций выхода нового фильтра, для соответствующих им плотностей вероятности найдена цепочка из интегродифференциального уравнения и формул тактовых пересчетов его решения, описан способ вычисления структурных функций фильтра методом Монте-Карло. Также построено гауссовское приближение к предлагаемому фильтру и его линеаризованное упрощение. Проведен сравнительный анализ этих приближений с аналогичными для классического фильтра. Точность и скорость оценивания линеаризованного фильтра проиллюстрирована двумерным примером, который подтвердил теоретические выводы.