CC BY
71
Об одном методе решения некорректно поставленных задач Драница Ю.П.
Судоводительский факультет МГТУ, кафедра высшей математики и программного обеспечения ЭВМ
Аннотация. Рассмотрен метод решения некорректно поставленных задач, основанный на теории калмановской фильтрации. Сделано обобщение этого метода на случай марковских процессов произвольного порядка и действия в системе коррелированного шума. Рассмотрен случай, когда для синтеза фильтра не имеется необходимой априорной информации, и предложены способы оценки его параметров, основанные на статистической обработке натурной информации. Разработан алгоритм параметрической адаптации фильтра к обрабатываемым данным. Метод предназначен для анализа и синтеза сложных технических и естественных систем.
Abstract. The method of the incorrectly stated tasks solution based on the Kalman theory filtration has been considered. The generalization of this method for the case of Markov processes of an arbitrary order in the correlated noise system has been made. The case, when for synthesis of the filter there is no necessary a priori information, has been considered, and the ways of the estimation of its parameters based on the statistical processing of the "on-location" information have been offered. The algorithm of parametrical adaptation of the filter to the processable data has been developed. The method is intended for the analysis and synthesis of complex technical and natural systems.
1. Общие положения
Во многих задачах интерпретации геофизических данных, зависимость fx), интересующая исследователя, не может быть явно измерена ни в одной точке пространства x. В то же время имеется модель, позволяющая связать искомую функцию с непосредственными измерениями другой функции F(x). Эта связь может быть отображена следующим операторным уравнением:
Afx, a) = F(x) + n, (1)
где: a - вектор параметров модели; А - оператор отображения; n - BeKTop ошибок наблюдений и моделирования.
Формально постановка задачи решения уравнения (1) состоит в следующем (Алгоритмы и программы..., 1984): задан непрерывный оператор, который взаимно однозначно отображает элементы f(x, a) метрического пространства Е1 в элементы F(x) метрического пространства Е2. Требуется найти решение операторного уравнения (1) в классе функций fx, a), если функция F(x) неизвестна, но имеются измерения y1,y2,...,yi функции F(x) в I точках x1,x2,...,xl.
Отметим, что "измерения" в контексте данной проблемы трактуются расширительно и могут включать не только какие-либо измерения физических величин, но и, например, результаты вычислительного эксперимента, решения прямых задач, некоторые эмпирические зависимости, графики, законы физики и любую другую априорную информацию, которую предполагается использовать при построении зависимостей. Оператор А уравнения (1) может представлять любое правило, по которому осуществляется отображение функциональных зависимостей. Мы же будем рассматривать отображения, осуществляемые посредством дифференциальных или интегро-дифференциальных уравнений.
Поставленную выше проблему часто называют задачей интерпретации результатов косвенных экспериментов. Постановка этой задачи может оказаться некорректной. Рассмотрим последнее более подробно.
Говорят, что решение операторного уравнения (1) устойчиво, если малая вариация правой части F(x) приводит к малому изменению решения. Это означает, что для всякого достаточно малого s найдется такое ад), что как только выполняется условие
P1(F(x),F(x+s)) < ад),
окажется выполненным и неравенство
P2(f(x)J(x+s)) <
Здесь Р1, Р2 означают расстояния, определяемые в метриках Е1 и Е2 соответственно. Говорят также, что задача решения операторного уравнения (1) поставлена корректно по Адамару, если его решение существует, единственно и устойчиво. Если не выполняется хотя бы одно из перечисленных требований, то говорят, что задача решения операторного уравнения поставлена некорректно.
Типичной некорректно поставленной задачей является, например, уравнение свертки:
Уравнение (2) является интегральным уравнением Фредгольма первого рода, отображающим множество непрерывных на отрезке [0,-т] функций f(t) на множество непрерывных на этом же отрезке функций S(t). Как показано в (Алгоритмы и программы..., 1984), решение уравнения (2) является некорректно поставленным в том смысле, что сколь угодно малым изменениям наблюденной функции S(t) могут соответствовать сколь угодно большие изменения зависимости f(t). Природа этого факта заключается в том, что интеграл в правой части уравнения (2) является сглаживающим оператором. Это сглаживание приводит к тому, что большие по амплитуде возмущения функции f(t), но сосредоточенные на малом участке пространства Е1, переводятся в малые возмущения функции S(t) (тем меньшие, чем меньше этот участок). При решении обратной задачи, чтобы объяснить малые возмущения функции S(t), приходится предполагать значительные отклонения функции f(t).
При решении задачи (2) на ЭВМ обычно прибегают к ее дискретизации. В этом случае дискретный аналог уравнения (2) представляется следующим выражением:
При сведении задач к моделям (2)-(3) обычно ситуация сильно идеализируется. В действительности, во-первых, наблюдения всегда осуществляются с некоторой погрешностью, и, таким образом, мы не располагаем точными значениями наблюдаемой зависимости. Во-вторых, наблюдения проводятся лишь на конечном множестве значений аргумента, и к тому же на дискретных его значениях (для модели (3)). Все это приводит к тому, что решение задач (2)-(3) не будет единственным, и различия между разными решениями будут большими. Другими словами, сильно отличающиеся решения ДО будут одинаково хорошо соответствовать наблюденной зависимости Е((). Следовательно, экспериментальные точки можно аппроксимировать множеством функций, по точности не уступающим истинной зависимости.
В современной практике имеется два подхода к решению этой проблемы. Первый метод заключается в разложении функций по заранее фиксированному конечному базису. В этом случае решение зависит от конечного числа параметров, и если число отсчетов превышает это число, то решение находится методом наименьших квадратов (МНК). Во втором случае рассматривается класс равномерно ограниченных функций. Тогда решение может быть найдено введением регуляризирующего функционала (Тихонов и др., 1990).
Отметим, однако, что рассматриваемые методы, хотя и позволяют получить устойчивые и однозначные решения, не гарантируют нахождения истинных функциональных зависимостей. Речь в этих случаях может идти только о построении эффективных моделей, т.е. нахождении зависимостей, которые достаточно точно аппроксимируют экспериментальные данные, но их интерпретация может не соответствовать реальной действительности. Дело в том, окончательное решение задачи обычно находится минимизацией того или иного функционала, который в общем случае многоэкстремален. Однако в настоящее время не существует методов, гарантирующих обнаружение глобального экстремума минимизируемой функции. Кроме того, постановка задач накладывает на изучаемые явления ряд достаточно серьезных ограничений (например, требование стационарности, нормальности действующих в системе шумов и т.п.), которые на практике выполняются довольно редко.
2. Основные принципы фильтрации Калмана
В связи с обозначенными выше проблемами, в последнее время все большее распространение получили методы решения некорректно поставленных задач, основанные на использовании априорной информации и позволяющие получить устойчивые и физически обоснованные решения.
S(t) = \f(t-T)u(T)dT.
(2)
S(k) = Ti f(k-i) U(i), i = 1,2,..., t.
(3)
Эти новые подходы имеют ряд принципиально важных отличий от рассмотренных выше. Классические методы, по сути дела, решают возникающие проблемы чисто формальным путем. Можно сказать, что они в большинстве случаев рассматривают изучаемый объект с позиций принципа "черного ящика", о котором известны только входные и выходные сигналы. При этом обычно предполагается стационарность исследуемого процесса. Основная возникающая при этом проблема заключается в определении импульсной характеристики исследуемой системы. Такой подход дает неплохие результаты при изучении простых систем, но в целом классические подходы к моделированию и анализу сложных технических или природных систем исчерпали себя.
Это проявилось еще в начале 50-х годов, и в результате предпринятых усилий возник принципиально новый подход к анализу и синтезу сложных систем, который часто называют методом, основанным на понятии пространства состояний. Теория этого метода заложены в работах Калмана -для дискретных систем и Калмана - Бьюси - для непрерывных (Браммер, Зиффлинг, 1982). Эти методы основаны на решении линейных дифференциальных уравнений и, следовательно, предполагают известной некоторую априорную информацию об изучаемом объекте. Они являются оптимальными по винеровскому критерию и позволяют снять ряд ограничений, присущих классическим подходам (например, не требуется предположения о стационарности). В настоящее время калмановский подход получил широкое распространение при решении различных навигационных задач (Hongzhuan, Hongyue, 1998), проблем автоматического регулирования и в других областях техники. Однако эта теория и ее возможные приложения, по нашему мнению, созвучны также широкому кругу геофизических задач. Кальмановская фильтрация имеет физически обоснованную теорию и законченный математический аппарат. Все это позволяет надеяться на получение принципиально новых свойств решений для такого важного класса проблем, как, например, некорректно поставленные задачи.
Динамические системы обычно описываются с помощью разностных, обыкновенных дифференциальных уравнений, или уравнений в частных производных, причем эти уравнения по своей природе могут быть как детерминированными, так и стохастическими. Отметим, что динамические факторы в этих уравнениях фигурируют в виде производных или разностных выражений. Известно, что любое обыкновенное дифференциальное или разностное уравнение n-ro порядка можно преобразовать в систему, состоящую из n дифференциальных или разностных уравнений первого порядка. Для этого вводят n так называемых фазовых переменных (или переменных состояния). Но в более общем случае переменные состояния определяются как множество таких переменных, которые вместе с известными входными сигналами необходимы и достаточны для нахождения отклика системы на внешнее воздействие. Из приведенного определения видно, что множество фазовых переменных не обладает единственностью. Число этих переменных определяет порядок системы, но наборы таких параметров для одной и той же модели могут быть различными. Физически этот феномен можно объяснить следующим образом. Поведение системы определяется множеством взаимосвязанных факторов, поэтому возможен выбор различных групп переменных, несущих идентичную информацию о поведении системы.
Задача моделирования системы в этом случае заключается в определении множества переменных состояния и их значений на заданном временном отрезке. Рассмотрим решение этой задачи для дискретного случая (Губанов, 1997). Это не умаляет общности рассуждений, так как наблюдения даже за непрерывными динамическими системами ведутся, как правило, дискретно, кроме того, имеется хорошо разработанный математический аппарат, обеспечивающий трансформацию дифференциальных уравнений в разностные. Будем считать, что наблюдения ведутся в равноотстоящие моменты времени to,tj,...,tk,tk+j,... В этом случае последовательность моментов наблюдения удобно обозначать их порядковыми номерами: 0,1,2,...,k,k+1,...
Пусть динамическая система описывается следующим разностным уравнением:
X(k+1) = A(k) X(k) + W(k), k = 0,1,... (4)
где: k - дискретный параметр времени; X(k) - неизвестный вектор состояния системы (его часто называют фазовым вектором) размерности n; W(k) - неизвестный n-мерный возмущающий вектор случайного входного воздействия; A(k) - детерминированная матрица размерности nxn (переходная матрица), элементы которой могут представлять собой функцию времени.
Соотношение (4) является разностным линейным уравнением первого порядка, которое связывает фазовые переменные системы для двух последовательных моментов времени и описывает марковский процесс первого порядка. Уравнение (4) называют дискретной моделью состояния системы. В частном случае, когда элементы матрицы А постоянны, модель (4) описывает стационарную
динамическую систему. Относительно вектора состояния обычно предполагается, что весь он или его часть могут быть недоступны непосредственному наблюдению. Матрица А(к) количественно описывает связь между последовательными состояниями системы. В общем случае она зависит от времени и должна быть получена на основе теоретических, экспериментальных или каких-либо других априорных соображений. Возмущающий вектор Щ(к), с одной стороны, характеризует шум, действующий на систему, а с другой - отображает неадекватность модели (4) описываемому процессу. Относительно вектора Щ(к) предполагается, что он центрирован, некоррелируем с фазовым вектором и представляет возмущения типа белого шума, для которого выполняются следующие соотношения:
Е{ Щ(к)} = 0, Е{ Щ{к)Х'(1)} = 0, (5)
Е{Щк)Щ(1)} = 2(к)4, (к, I = 0,1,2,...), (6)
где: - символ Кронекера (^ = 0, если к Ф I и 5к,- = 1, если к = /); Е{} - оператор вычисления математического ожидания; значок "' " обозначает операцию транспонирования.
Будем также предполагать, что динамическая система (4) наблюдаема, т.е. имеет пг-мерный выходной сигнал У(к), который может быть измерен в любой момент времени и который несет в себе информацию о состоянии системы Х(к). Будем считать, что между наблюденным выходным сигналом и вектором состояния имеется следующая линейная связь:
У(к) = С(к) Х(к) + У(к), (7)
где: С(к) - заданная матрица размерности пхпг; У(к) - неизвестный вектор погрешности измерений (невязок) размерности пг.
Выражение (7) называется моделью наблюдений. Отметим, что вектор погрешностей У(к) в общем случае содержит как ошибки измерительного тракта, так и невязки, связанные с неадекватностью модели (7) измеряемым данным. Относительно него делаются предположения о центрированности, некоррелируемости с фазовым вектором и вектором возмущений Щ(к):
Е{У(к)} = 0, Е{У(к)Х'(к)} = 0, Е{У(к)№(к)} = 0. (8)
Предполагается также, что он представляет реализацию типа белого шума и, следовательно, имеет ковариационную матрицу следующего вида:
Е{У(к)У(Г)} = Я(к)4. (9)
Для дискретной системы обычно предполагается, что возмущение Уф представляет собой ступенчатую функцию типа телеграфного сигнала, постоянную на каждом интервале дискретизации [4,4+Л длиной Т, где Т - шаг дискретизации.
Отметим, что выполнение требований (5)-(6), (8)-(9) для данной методики является принципиальным. Физически они означают, что возмущения и шумы, действующие в системе, являются белыми (некоррелируемыми). Однако эти требования не являются жесткими, так как их всегда можно обойти за счет некоторого усложнения модели. Действительно, пусть действующий в системе шум Щ(к) является цветным. Такой шум можно представить как выход некоторой дополнительной динамической системы, возбуждаемой заведомо некоррелируемым белым шумом и(к) (Сейдж, Меле, 1976):
Р(к) = М(к-1) Р(к-1) + и(к),
Щк) = Н(к) Р(к) + Т(к), (10)
где не только и(к), но и Т(к) - белые некоррелированные шумы.
Рассмотрим теперь обе динамические системы (4), (7) и (10) совместно, т.е. введем новый фазовый вектор
Х(к)
Дк) =
Р(к)
а также следующие матрицы:
т =
Л(к) Н(к) 0 М(к)
Цк) =
Е 0 0 Е
где Е - единичная матрица.
Легко убедиться, что вектор Z(k) удовлетворяет уравнению
Z(k) = Т(к-1) Z(k-1) + Ц(к) Щ(к),
(12)
где
Щ(к) =
Т(к)
и(к)
(13)
представляет собой составной вектор белых шумов.
Полученное уравнение (12) полностью совпадает по форме с уравнением (4), поэтому его можно назвать приведением модели динамической системы с цветным шумом Щ(к) к модели с белым шумом Щ(к). Благодаря процедуре приведения снимается условие о том, что возмущение динамической системы должно быть белым шумом, и вся теория калмановской фильтрации, построенная при этом условии, остается в силе за счет некоторого усложнения модели этой системы (Губанов, 1997).
Классическая задача теории пространства состояний заключается в определении несмещенной (с математическим ожиданием, равным истинному значению оцениваемых параметров) оценки Хл(к) вектора Х(к) по заданным измерениям У(0),У(1),...,У(к). Требуется на основе имеющейся последовательности данных построить такую линейную несмещенную оценку фазового вектора Хл, чтобы дисперсия ошибок
Х°(к) = Х(к) - Хл(к),
т.е. диагональные элементы матрицы апостериорных автоковариации
Р(к) = Е{Х°(Х°)'}
были минимальными, при условии выполнения соотношений (5), (6), (8), (9). Именно такая оценка называется оптимальной.
Как и при решении дифференциальных уравнений, метод требует задания начальных условий, в качестве которых выступают: ХЛ(0) - оценка начального вектора состояния системы; Р(0) - матрица ковариаций ошибок оценки вектора начального состояния. Кроме того, предполагаются известными ковариационные матрицы Ц(к) возмущающего вектора Щ(к) и ковариационной матрицы Я(к) вектора ошибок измерения У(к). В этих предположениях, процесс эволюции во времени оптимальной по Калману оценки фазового вектора описывается следующими рекурентными соотношениями:
Х°(к) = Л(к-1) Х\к-1), (14)
Х(к) = Х°(к) + К(к) [У(к) - С(к) Х°(к)], (15)
где ^л(0) - задано; к = 1,2,...
Необходимые матрицы определяются также рекурентными соотношениями
Р(к) = Л'(к-1) Р(к-1) Л(к-1) + 2(к-1), К(к) = Р(к) С(к) [С'(к) Р(к) С(к) + Я(к)]-1, Р (к) = (Е - К(к)С'(к)] Р(к), Рл(0) = Р(0), (16)
где значком "-1" обозначены операция вычисления обратной матрицы; Е - единичная матрица соответствующей размерности.
Матрица К(к) называется матрицей усиления фильтра Калмана. Матрица Р(к) характеризует эволюцию ошибки оценки фазового вектора. Работа алгоритма заключается в следующем. По
соотношению (14) выполняется предварительная и приближенная экстраполяция фазового вектора на один шаг вперед. Затем экстраполированная оценка по (15) уточняется с учетом новых данных наблюдений. Полученная таким образом улучшенная оценка используется для прогноза вектора состояния на следующем шаге и т.д. Аналогичным образом по соотношениям (16) экстраполируются, а затем и уточняются ковариационная матрица ошибок оценки вектора состояния и матрица усиления фильтра. Как видно из соотношения (16), стабилизация решений достигается за счет учета характеристик действующих в системе шумов.
Рекуррентный алгоритм (14)-(16) называется фильтром Калмана. Он состоит из двух частей. Собственно фильтром является та часть алгоритма (формулы (14) и (15)), где производится экстраполяция и уточнение фазового вектора. Вторая часть алгоритма, которая называется синтезом оптимального фильтра, заключается в вычислении матрицы усиления фильтра и ковариационных матриц ошибок экстраполяции и оценивания (формулы (16)), которые не зависят от данных наблюдений. Поскольку матрицы А, С, 2, Я предполагаются известными, эта часть алгоритма может быть рассчитана заранее и запомнена.
3. Обобщение теории Калмана
Из приведенных выше соотношений следует, что при синтезе оптимального фильтра требуется значительный объем априорной информации. Прежде всего это начальное значение фазового вектора АА(0) и матрица ковариаций ошибок этой оценки Р(0). Принципиальным является знание ковариационных матриц действующих в системе шумов. Адекватное задание этой информации определяет не только точностные свойства фильтра, но и его устойчивость. Другими словами, некорректное определение этой информации может привести не только к большим ошибкам алгоритма, но и даже к расхождению оценок фазового вектора. Однако для большинства практических задач эта информация априори не известна, и требуется ее оценка. Для этого, как правило, прибегают к дополнительной обработке натурных данных. Рассмотрим этот вопрос более подробно.
Предположим, что не имеется никакой информации о начальном состоянии динамической системы в момент времени к = 0. В этом случае выполним, начиная с момента к = 0, число наблюдений /, превышающее число фазовых переменных п. Представим эти данные системой линейных уравнений с постоянными параметрами
и решим ее по МНК, полагая, что матрица ковариаций ошибок измерений известна априори. Причем эта матрица необязательно должна быть диагональной, т.е. шум У(к) на интервале времени [/0, может быть цветным. В результате получим некоторую усредненную на этом интервале оценку ХЛ(к) и апостериорную ковариацию Р(к) ее ошибок. Примем полученные оценки фазового вектора и матрицу ковариации ее ошибок в качестве начального состояния системы в момент времени к.
Теперь с помощью уравнения (14) можно построить приближенный прогноз фазового вектора на следующий момент наблюдений:
который не учитывает возмущения Щ(к). Оценка (18) получена благодаря знанию переходной матрицы А(к-1) с учетом предшествующих наблюдений. Таким образом, возможен старт метода и без знания начального состояния фазового вектора.
Рассмотрим теперь случай, когда из априорных соображений не удается построить переходную матрицу А(к) модели (4), но априорно известна ковариационная матрица фазового вектора б(к). В практических приложениях такая ситуация возникает довольно часто, когда фазовые переменные могут быть оценены или непосредственно измерены, например, по другому независимому материалу, решением прямой задачи. Если такие данные имеются, то корректная оценка ковариационной матрицы может быть выполнена, например, по методике, изложенной в работе (Драница, 1999). Рассмотрим теперь задачу построения дискретной динамической модели для коррелированного процесса Х(/) с заданной матрицей автоковариации б(к):
У(У) = С() Х + У(/), ] = 1,2,.../ (/ > я)
(17)
Х°(к) = А(к-1)Хл(к-1),
(18)
С(к2,к1) = Е{Х(к2)Х' (к1)}.
Если найденная матрица удовлетворяет следующим условиям:
а) для всех к матрица дисперсий О(к) = Е{Х(к)Х(к)} задана, конечна и положительно определена;
б) для всех к1 < к2 < к3 имеет место соотношение
О(к3,к2)О'1(к2)О(к2,к1) = О(к3,к1), (20) то согласно работе (Губанов, 1997), переходная матрица может быть представлена следующим выражением:
Л (к) = О(к-1,к)О(к)-1. (21)
Отметим, что условие (20) справедливо только для цепей Маркова, т.е. случайных последовательностей без последействия. Как показано в работе (Драница, 1999), соотношение (21) выполняется также и для одношагового линейного предсказателя, который в рассматриваемом контексте может интерпретироваться как марковский процесс первого порядка.
Найдем теперь матрицу ковариаций Q(k) вектора возмущений Щ(к), если известна последовательность ковариационных матриц фазового вектора Х(к) (19). Для этого рассмотрим матрицу дисперсий О(к+1). Подставим выражение (4) в уравнение (19), после чего получим:
О(к+1) = Е{Х(к+1)Х(к+1)} = Е{(Л(к)Х(к) + Щ(к))(Л(к)Х(к) + Щ(к))'}. (22)
Поскольку Х(к) и Щ(к) некоррелируемы, математическое ожидание их произведения равно нулю (см. выражение (5)), поэтому формула (21) упрощается до следующего выражения:
О(к+1) = Л(к)О(к)Л'(к) + О(к). (23)
Поскольку матрицы дисперсий О заданы, а переходная матрица А известна, поэтому непосредственно из (23) находим требуемую ковариационную матрицу возмущений:
Q(k) = О(к+1) - Л (к) О(к) Л'(к). (24)
Допустим теперь, что требуется оценить матрицу ковариаций Я(к) вектора ошибок измерения У(к). Поступим аналогично предыдущему случаю. Определим математическое ожидание левой и правой частей выражения (7).
Е{У(к) (У(к))'} = Е{(С(к) Х(к) + У(к)) (С(к)Х(к) + У(к))'}. (25)
Учитывая некоррелируемость векторов Х(к) и У(к), будем иметь
Е{У(к)У(к)} = С(к) Е{Х(к)(Х'(к)} С (к) + Е{У(к)У(к)} = С(к) О(к) С (к) + Я(к).
Из последнего уравнения непосредственно следует
Я(к) = Е{У(к)У'(к)} - С(к) О(к) С'(к). (26)
Таким образом, всю последовательность матриц Л(к) можно определить из соотношения (21), после чего по формулам (24), (26) определяются и последовательности матриц Ц(к), Я(к). Начальное значение фазового вектора Х(к) и ковариацию ошибки его оценки Р(к) определяем согласно алгоритму (17)-(18). Тем самым процесс синтеза оптимального фильтра оказывается полностью определенным. Полученная при этом модель описывает динамическую систему, способную синтезировать из белого шума Щ(к) Марковский коррелированный процесс Х(к) с заданными статистическими свойствами, т.е. ковариации (19), удовлетворяющие сформулированным выше условиям (а) и (б). При этом, однако, все требуемые для синтеза фильтра матрицы должны вычисляться в оперативном режиме. Преимуществом такой фильтрации является его адаптация к имеющимся данным. На начальном этапе работы алгоритма можно ограничиваться довольно грубыми оценками требуемой информации. По мере обработки свежих данных эти оценки могут уточняются, адаптируя таким образом модель к конкретному материалу.
Как указывалось выше, модель (4), (7) описывает марковский процесс первого порядка. Иначе говоря алгоритм "помнит" только предыдущее состояние системы и "не помнит" всей ее предыстории. Однако реальные сложные системы не обязательно являются такими, а значит, их динамическая модель может иметь порядок выше первого. Это, в свою очередь, означает, что для адекватного описания системы требуются значения не одного, а нескольких предыдущих состояний фазового вектора. Обобщим алгоритм Калмана и на этот случай.
Допустим, что динамическая система описывается линейным уравнением с произвольной глубиной учета ее предыстории, т.е. может быть представлена следующим разностным уравнением:
Х(к+1) = Т,А,(к)Х(к-,) + ^(к+1), , = 1,2,...,г; т> 1, (27)
где субматрицы А,(к) имеют размерность яхя,, = 1,2,...,ги следующую структуру:
а12, а1я,
А,(к) = а21, а22, а2я, (28)
ая1, ая2, аяя,
Отметим, что субматрица А, представляет вклад в текущий прогноз фазового вектора, сдвинутого на, временных тактов.
Модель (27) представляет, с одной стороны, уравнение многомерной линейной регрессии, а с другой - описывает марковский процесс г-ого порядка. Часто глубину учета предыстории процесса называют лагом модели. Чтобы можно было воспользоваться описанной выше теорией калмановской фильтрации, необходимо привести модель (27) к уравнению первого порядка. Для этого введем следующий составной фазовый вектор:
Х°(к) = Х(к) иХ(к-1) II...II Х(к-т),
(29)
где знаком "и" обозначена операция объединения множеств.
Тогда, с учетом преобразования (29), уравнение (27) может быть записано следующим образом:
Х°(к+1) = А (к) Х°(к) + ^°(к+1),
(30)
где Х°(к) - составной фазовый вектор имеет размерность яхт; А(к) - переходная матрица размерности (яхяхг)х(яхяхт); - обобщенный вектор шума размера ях т.
Введенный фазовый вектор Х° представляет просто упорядоченную последовательность (т-1) исходных фазовых векторов. Обозначим вектором Ак, к-ую строку субматрицы (27). Тогда составная матрица А(к), с учетом преобразования (29) и обозначений (28), будет иметь следующую структуру:
А (к)
А11 А12 ■■■ А1(т-1) А21 А22 ■■■ А2(т-1)
АА
Ан
А2т
Аш(т-1) Атг
*т1 ■**т2
Е 0 ... 0 0 0 Е ... 0 0
0 0 ... Е
0
(31)
где Е, 0 - соответственно единичная и нулевая матрицы размерности яхя.
С учетом временной некоррелируемости вектора ошибок №(к) с исходными фазовыми векторами, обобщенный вектор №°(к) может быть выражен следующим образом:
= У О,
(32)
где О - нулевой вектор размера ях(г-1).
В этом случае ковариационная матрица вектора (2°) очевидным образом связана с соответствующей матрицей уравнения Маркова первого порядка (6):
б°(к) =
0,к) о ... о о о ... о
о о ... о
где о - нулевые матрицы размерности пхп.
И, наконец, с учетом преобразования (29), модель измерения (7) трансформируется к следующему виду:
У(к) = С°(к) Х°(к) + У (к) с блочной матрицей С°(к) размера пгхп(тхп) и следующего вида:
(34)
С°(к) =
С(к) О ... О О О ... О
О О ... О
(35)
где О - нулевая матрица размера пгхп.
Ясно, что ковариационная матрица вектора ошибок наблюдений У(к) не изменится, т.е. Я° = Я. Полученные уравнения (30), (34) полностью совпадают по форме с уравнениями (4), (7), поэтому этот алгоритм можно назвать приведением марковского процесса г-ого порядка к марковскому процессу первого порядка. Благодаря этой процедуре снимается условие, что описываемая система должна быть марковской первого порядка, и вся теория фильтрации Калмана, построенная при этом условии, остается в силе за счет некоторого усложнения модели системы. Необходимые для этого матричные преобразования приводятся в формулах (31), (33), (35).
Отметим далее, что представление шумовой компоненты соотношением (33) не является единственно возможным. В более общем случае она может иметь блочно-диагональный вид с ненулевыми матрицами на главной диагонали. Последнее означает, что на входе системы фактически может действовать коррелируемый по времени цветной шум, что значительно расширяет возможности моделирования.
4. Пример построения фильтра Калмана
Проиллюстрируем теперь метод на примере решения уравнения (3). Для простоты примем, что входом и выходом этого процесса описываются одномерными временными рядами. В этом случае соотношение (3) описывает некоторую одномерную систему с известной переходной функцией П(к) с числом отсчетов т. Измеряемым выходом системы служит скалярная величина S(k), а входом является коррелированный скалярный шум /(к), который требуется оценить по сделанным измерениям. В этих предположениях соотношение (3) будет представлять модель одномерной системы, которая описывается марковским процессом г-ого порядка. Для того, чтобы можно было воспользоваться теорией Калмана, преобразуем уравнение (3) в марковскую систему первого порядка. Определим следующий фазовый вектор размерности -л
Х°(к) = /(к-1) и/к-2) и...иАк-т\ (36)
Нам необходимо связать два последовательных состояния фазового вектора Х°. Если об этом имеется какая-либо априорная информация, примем ее к рассмотрению, т.е. тем или иным образом определим матрицу А соотношения (38). Будем, однако, предполагать, что такая информация отсутствует. Принимая во внимание, что входной вектор представляет реализацию типа белого шума, определить связи двух последовательных состояний можно, по крайней мере, на основе двух допущений. В первом случае будем полагать, что матрица А является единичной, т.е. будем считать, что Х°(к) = Х°(к-1). Последнее означает, что входное воздействие представляет некоторый стационарный процесс. Во втором случае будем считать, что текущее значение фазового вектора является средним взвешенным его предыдущих состояний. В такой постановке требование стационарности входного воздействия можно не предполагать, что является более общим случаем по сравнению с первым
вариантом. Поэтому рассмотрим второй подход более подробно. В этом случае переходная матрица динамической системы будет иметь следующий вид:
А =
Ь Ь2 ■ ЬТ-1 Ьт 1 0 ... 0 0
0 0 ... 1 0
(37)
где Ь1 = Ь2 = ... = Ьт= 1/(г-1) - весовые коэффициенты.
Матрица (37) имеет простую структуру, в которой первой строкой является строка весовых коэффициентов, под главной диагональю стоят единицы, а все остальные элементы имеют нулевое значение. При сделанных предположениях уравнение, описывающее эволюцию фазового вектора, будет иметь следующий вид:
Х°(к) = А Х°(к-1) (38)
с измерительным трактом следующего вида:
= Т, Х°(к-1)П(1). (39)
Перейдем теперь к оценке начального значения фазового вектора. Если имеется какая-либо априорная информация о начальном состоянии фазового вектора, примем ее к сведению, иначе выполним необходимое число измерений £(/), составим линейную систему (17) и, решив ее по МНК, определим значение этого вектора.
Рассмотрим теперь решение этой же проблемы, но с частотных позиций, сначала для непрерывного времени. В этом контексте свертку (2) можно рассматривать как общее представление произвольного линейного стационарного оператора с конечной временной импульсной характеристикой £/(/). Применим преобразование Фурье к обеим частям равенства (2). Известно, что преобразование Фурье от свертки функций равно произведению спектров этих функций, в этом случае получим
ЗД = И(ф)/», или (40)
/(ю) = ЗД/И(Й>), (40')
где/(&), £(®), Н(а>) - спектры соответственно входного, выходного сигналов и импульсной функции и.
Выполним обратное преобразование Фурье для соотношения (40'), т.е. вновь перейдем во временную область:
М = (1/2я) I (ЗД / И(®)) ехр(1й*) (41)
Полученное решение (41) имеет особенность при И(а>) = 0, т.е. при нулевом значении знаменателя. В этих точках может нарушаться единственность и однозначность решений. Для устранения таких дефектов можно использовать регуляризацию по норме (Тихонов и др., 1990). Для этого при интегрировании (41), будем использовать следующую частотную характеристику
/» = ЗД И\а) / (|И(®)|2 +02), (42)
где И - число, комплексно-сопряженное к И; а > 0 - параметр регуляризации; |И| - модуль И.
Параметр регуляризации находится методами оптимизации или определяется из априорных соображений. Отметим, что, согласно теории, регуляризирующая добавка в (42) представляет дисперсию белого шума, воздействующего на систему. Уравнение (42) всегда имеет значение знаменателя, отличное от нуля, и, следовательно, однозначно разрешимо (хотя и с потерей точности). Обозначим получаемые тем или иным методом решения символом Х*.
Переход к дискретному времени не вносит в этот алгоритм принципиальных затруднений, поэтому и для этого случая проблема решается аналогичным образом. Добавим только, что полученные на этом этапе оценки решения Х* могут быть использованы при идентификации модели (4) и уточнения параметров переходной матрицы А.
Таким образом, для данной задачи имеется несколько способов определения начального состояния фазового вектора. В тоже время спектральный подход имеет некоторые преимущества по сравнению с временным. Например, из чисто формальных соображений неясно, имеет ли уравнение (2) решение, а если имеет, то сколько. Переход в спектральную область сводит решение проблемы к ряду операций, всегда имеющих однозначный выход (за исключения случая, когда И(а>) = 0, но и в этом случае имеется способ получить однозначное решение). Следовательно, появляется возможность не только вычислить явное решение, но и исследовать области его существования и единственности. Кроме того, прямое и обратное спектральные преобразования могут быть выполнены алгоритмом быстрого преобразования Фурье, что резко снижает вычислительные затраты и является необходимым при массовой обработке материала.
Далее необходимо оценить ковариационные матрицы действующих в системе шумов. Отметим, что матрица шумов измерительного канала (39) трансформируется в скаляр, ковариация которого при известной матрице 2 может быть непосредственно оценена по приведенной выше методике (алгоритм (26)). Рассмотрим теперь оценку матрицы 2. Если временной ряд /(к) может быть оценен по каким-либо другим, независимым наблюдениям, то оценку шумов модели (38) также можно сделать непосредственно 2 = Е{Х*/}. Будем, однако предполагать, что какой-либо дополнительной информации по шумам в системе не имеется. В этом случае требуемую оценку можно выполнить, по крайней мере, двумя способами.
Зададим какие-либо произвольные значения шумов, действующих в системе и рассмотрим квадрат невязки
я = Е/ т - SЛ(/))2. (43)
Будем с помощью какой-либо оптимизационной процедуры подбирать параметры ковариационной матрицы шума 2(к) для достижения минимума функционала (43). В настоящее время процедуры оптимизации входят во многие математические и статистические пакеты прикладных программ, т.е. оптимизация является достаточно рутинной операцией, поэтому не будем рассматривать конкретную реализацию решения этой проблемы. Параметры полученных таким образом матриц примем за действующий в системе шум. Недостатками такого подхода являются значительные вычислительные издержки, требующиеся в процессе оптимизации параметров шума.
Рассмотрим теперь другой способ оценки параметров шума с меньшими вычислительными издержками. Зададим начальные оценки параметров шума следующим образом:
2 = с?Е, (44)
где Е - единичная матрица размера гх т, а2 - дисперсия фазового вектора /(к).
Требуемая дисперсия, входящая в формулу (44), может быть получена, например, по ранее полученным оценкам фазового вектора Х*. Отметим, что с позиций МНК оценка (44) имеет место для ортонормированного базиса и единичной порождающей матрицей Е. Эта оценка оптимальна, когда параметры линейной модели постоянны, а их невязки представляют собой случайные ошибки наблюдений, которые центрированы, нормально распределены, независимы и равноточны. Можно принять эту оценку за начальную, а затем выполнить описанную выше процедуру оптимизации. В этом случае простейшая матрица ковариаций (44) будет заменена на более общую диагональную матрицу 2 = diag(a;2), составленную из различных дисперсий отдельных измерений. Оценка аналогичного вида может быть получена, если рассчитать ковариацию фазового вектора Х*. В линейных пространствах таким матрицам соответствует ортогональный, но не нормированный базис. Другой способ оценки матрицы 2 заключается в вычислении ковариации 2 = Е{(ХЛ-Х*)(Х"'-Х*)'}. Принимая во внимание, что при оценке вектора Х*, шум модели явно не учитывается, можно предполагать, что такой способ оценки более адекватен физике явления.
Определим ковариацию вектора ошибки оценки следующим образом
Р = 8Е, (45)
где Е - по-прежнему единичная матрица; 3 - скалярная величина, представляющая собой верхнюю границу значений фазового вектора, т.е. 5= 3 а. Если исходить из "правила трех сигм", то оценка (45) предполагает, что никакой априорной информации о точности начального оценивания не имеется. Как и
для других оценок, она может быть улучшена соответствующей процедурой оптимизации. Очевидно, что рассмотренный пример синтеза модели одномерной системы без труда может быть обобщен и на многомерный случай.
Отметим следующие важные свойства калмановской фильтрации. Рассмотрим последовательность оценок фильтрации Хл(0), Хл(1), Хл(2), ..., XA(k) стохастического параметра динамической системы. Эта последовательность обладает кумулятивными свойствами в том смысле, что их точность неоднородна и растет с ростом k, дисперсия ошибок оценивания, соответственно, уменьшается, а в пределе они становятся константами. Если кумулятивные свойства получаемых решений нежелательны, рекомендуется их сгладить способом обратной фильтрации (Губанов, 1997).
Решение уравнения свертки (3) методом фильтрации Калмана было запрограммировано и протестировано на натурных и синтетических данных. Экспериментами на этой модели установлено, что достаточно высокая стабилизация решения достигается уже на временах t «г. Из этого следует, что если выполнить планирование эксперимента таким образом, чтобы можно было пренебречь этим начальным участком решения, то требования к точности оценки начального значения фазового вектора и дисперсии этого оценивания становятся непринципиальными.
5. Заключение
Таким образом, разработанная методика позволяет моделировать нестационарные многомерные линейные динамические системы и процессы произвольного порядка, используя теорию фильтрации Калмана. Предлагаются способы оценки параметров синтезируемого алгоритма и характеристик действующих в системе шумов, основанные на статистической обработке натурных данных. Имеется возможность описания систем как с белым, так и цветным шумами на входе. Важным достоинством разработанного метода является возможность его адаптации к натурным данным.
Метод предназначен для анализа и синтеза сложных технических и естественных систем. Известно, что калмановская фильтрация обладает рядом полезных свойств: хорошая помехозащищенность и устойчивость, а также оптимальность по винеровскому критерию. Расширение возможностей этого метода должно улучшить качество моделирования сложных процессов и систем. Следует ожидать, что разработанная методика может оказаться особенно полезной при решении обратных задач.
Литература
Hongzhuan Qiu and Hongyue Zhang. Federated filter and its application to integrated navigation system.
Second International Symposium on Inertial Technology in Beijing, October 27-29, p.18-23, 1998.
Алгоритмы и программы восстановления зависимостей. Под ред. Вапника В.Н. М., Наука, 816 е., 1984.
Браммер К., Зиффлинг Г. Фильтр Калмана - Бьюси. М., Наука, 552 е., 1982.
Губанов B.C. Обобщенный метод наименьших квадратов (теория и применение в астронометрии). СПб, Наука, 319 е., 1997.
Драница Ю.П. Об одном методе моделирования нестационарных динамических систем и процессов.
Вестник МГТУ, т.3, № 1, с.55, 2000. Сейдж Э., Меле Дж. Теория оценивания и ее применение в связи и управлении. М., Связь, с.53, 1976. Тихонов А.Н., Гончарский A.B., Степанов В.В., Ягола А.Г. Численные методы решения некорректных задач. М., Наука, 232 е., 1990.
