Анализ вычислительных трудностей при реализации фильтра Калмана для оценки полного электронного содержания ионосферы Текст научной статьи по специальности «Математика»

22

УДК 550.388.2

А. Ф. Лаговский, А. В. Мышерин, И. И. Шагимуратов

АНАЛИЗ ВЫЧИСЛИТЕЛЬНЫХ ТРУДНОСТЕЙ ПРИ РЕАЛИЗАЦИИ ФИЛЬТРА КАЛМАНА ДЛЯ ОЦЕНКИ ПОЛНОГО ЭЛЕКТРОННОГО СОДЕРЖАНИЯ ИОНОСФЕРЫ

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

The article examines Kalman filter implementation for estimation of the total ionospheric electron content. The major methods allowing to avoid an uncontrollable spread of mistakes are examined. Basic approaches permitting to get numerical stability at construction of the filter are considered. The most effective numerical methods are identified for this model.

Диагностике состояния ионосферы уделяется большое значение в различных программах геофизических и космических исследований. Свойства ионосферы таковы, что она изменяет траекторию радиосигналов от спутников, являясь источником ошибок измерений, проводимых на Земле. Измерение и анализ этих ошибок позволяет оценивать некоторые параметры ионосферы, которые помогают больше понять её природу, построить наиболее точную её модель.

Задача оценки полного электронного содержания (ПЭС) ионосферы уже рассматривалась нами [1], поэтому приведем лишь основные формулы модели ПЭС в форме Георгиади:

h = 12 (T -14),

Nh(t) = a0 + a1 cosh + a2 cos2h +... + an cosnh + (1)

+ b1 sin h + b2 sin2h +... + bn sin nh + с1Аф + c2 A<ph + c 3Аф2.

Здесь LT — местное время, Аф — разности между широтой станции и подионосферной точки.

В нашем случае n равно шести. Для каждого момента времени уравнение (1) эквивалентно системе уравнений, в которой неизвестными являются 16 коэффициентов разложения ионосферной задержки, аппаратурные поправки АП,С считаем известными.

Рассмотрим общую модель фильтра Калмана [2]

Вестник РГУ им. И. Канта. 2007. Вып. 10. Физико-математические науки. С. 22 — 27.

Динамическая модель системы (уравнение изменения параметров модети) хк = фк_1хк_1 + Щ-1.

Модель наблюдаемых значений хк = Нкхк + ук, где т ~ N(0, Qk) и Юк~ N(0, Як) — аддитивные помехи.

Нача.льные условия М(.го) = *о, М(*о*от) = Ро.

Условия независимости М() = 0 .

Д.ля .любых к к ] Хд(-)=Фа_1 хк_г (+).

Априорная оценка матрицы ковариации ошибок:

Рк (-) = ф к-1 Рк-1(+)ф!-1 + Qk-l.

Апостериорная оценка вектора параметров модели:

Хк(+)=%к(-) + Кк[Хк — -НЛ(-)].

Апостериорная оценка матрицы ковариации ошибок:

Рк (+) = [1 - КкНк ]Рк (-).

Коэффициент усиления (ядро) фильтра Калмана:

Кк = Рк (-Н [ НкРк (-Н + Як ]-1.

Итак, вычислительный алгоритм должен иметь следующий виц.

1. Вычисляем Рк(-), используя Рк - 1(+), Фк - 1 и Qk - 1.

2. Вычисляем Кк, используя Рк(-) (из пункта 1), Нк и Як.

3. Вычисляем Рк(+), используя К (из пункта 2), и Рк(-) (взятого из пункта 1).

4. Рекурренгно вычисляем величины Т/ (+), используя значение величины Кк (из пункта 3), заданную начальную оценку %о и входное значение Хк.

В нашем случае Фк — единичная матрица (так как параметры модели можно считать постоянными, по крайней мере, в течение суток), а строки матрицы Нк формируются из уравнения (1).

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

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

Главная причина численных трудностей в реализации обычного фильтра Калмана как с точки зрения объема вычислений, так и с точки зрения вычислительных ошибок — решение матричного уравнения Риккати.

23

24

Уравнение имеет виц:

ф л Ф1 + Ок = Рк+1-

Доказано, что если матрица перехода состояния Фк несингулярна и

Рк (-) = АкБ-1

есть несингулярное матричное решение дискретного уравнения Рик-кати для момента времени Ьк, то

Рк+1(-) = Ак+1Б-+1 есть решение в момент времени 1к+ъ где

Ак+1 & I" ф-Т 0 1

_ Бк+1 _ I 0 0 фк _

Щя-1 Нк

I

Ак

Бк

ф к + Ок Фк1'НкЯ-1 Нк Ок Ф-Т' 1 к 1

Ф кТНТкЯ-1Нк ФкТ _ 1 к |

Неконтролируемое распространение ошибок в решении уравнения Риккати — главная причина деградации в работе фильтра. В нашем случае ядро фильтра очень быстро растет, что приводит к переполнению, и фильтр не успевает обработать даже 24-часовой интервал.

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

Наиболее стабильные численные реализации фильтра Калмана используют одну или более из перечисленных ниже техник [2; 6; 7].

1. Разложение на множители матрицы ковариации неопределенности оценки состояния Р (разложение Холецкого).

2. Разложение на множители матрицы ковариации шума измерения Я (эти методы фактически «декоррелируют» компоненты вектора шума).

3. Взятие симметричных матричных квадратных корней элементарных матриц. Симметричная элементарная матрица имеет форму

I - оии1, где I — единичная матрица п х п, а — скаляр и V — вектор. Симметричный квадратный корень из элементарной матрицы — также элементарная матрица с тем же самым V, но другим значением а.

4. Основные разложения матриц на множители как произведения треугольной и ортогональной матриц.

Методы ОЯ-разложения были первоначально развиты для нахождения устойчивых решений систем линейных уравнений. Матричные множители в произведении ортогональной матрицы О и треугольной матрицы Я. В применении к фильтрации Калмана необходим только треугольный множитель. О и Я уже имеют специальные значения в фильтрации Калмана. Существуют два основных метода построения ОЯ разложения:

— вращения Гивенса, одновременно оперируют одним элементом;

— отражения Хаусхолдера, одновременно оперируют одним столбцом.

В дополнение следует отметить ортогонализацию Грамма — Шмидта как другой общий метод для разложения на множители матрицы общего вида. Он позволяет представить матрицу в виде произведения прямоугольной матрицы с ортонормированными столбцами и треугольной матрицы. Для фильтра Калмана важен только треугольный множитель.

5. Алгоритмы одноранговой модификации. Одноранговая модификация симметричной положительно определенной п х п матрицы М имеет вид М. ± vvT, где V — п-мерный вектор.

Рассмотрим подробней те методы, которые возможно применить для улучшения численных свойств фильтра Калмана нашей модели. Если М — положительно определенная симметричная матрица, то разложение Холецкого для нее выглядит следующим образом:

ССТ=М.

Здесь С — нижняя треугольная матрица. Такое разложение определяется неоднозначно и имеет вычислительную сложность порядка п3 /3 (п — порядок матрицы).

Также для симметричной положительно определенной матрицы можно построить модифицированное разложение Холецкого:

М=Шит.

Здесь и — верхняя треугольная, а О — диагональная матрица. Сложность алгоритма п3 /6.

Рассмотрим, как с помощью таких разложений можно произвести декорреляцию ошибок наблюдения. Ошибки появляются из уравнения

2 = Нх + V ,

где М(vvT) = Я — недиагональная матрица. Тогда скалярные компоненты 2 могут обрабатываться последовательно как скалярные наблюдения со статистически независимыми ошибками наблюдения. Впрочем, матрица Я всегда может быть представлена в виде

я = иои,

где деШ = 1, и, следовательно, она будет всегда иметь обратную матрицу. Причем необязательно вычислять Ш-1 для декорреляции, достаточно переопределить наблюдения как

2 = и^2 = и-1(Нх + V) = (и^ Н)х + и^ (V) = НX + V,

Я' = М( V' V'Т) = М(и-1 тТ (иТ )-1) = и-1 М( шТ )иТ-1 =

= и -1(иоиТ )иТ-1 = О .

Таким образом, для декорреляции ошибок наблюдений необходимо решить систему уравнений

и2 = 2,

иН ' = Н.

25

26

При этом не нужно обращать и, так как обе системы имеют верхнюю треугольную матрицу.

Процедура декорреляции в несколько другом виде (с помощью ор-тогонализации Грамма — Шмидта) показана в [5]. Для нашей задачи приведенный выше метод выглядит предпочтительнее.

Для расчета коэффициента усиления фильтра Калмана необходимо обращать матрицу ИРИТ+Я. Для нашей модели эта матрица является жесткой, поэтому такое обращение вызывает наиболее существенные вычислительные ошибки. Избежать этого можно следующим образом: представим ИРИТ+Я = иои1. Тогда

[иоиТ ] • [ИРИТ + я]-1 И = И .

То есть имеем систему

ИБЦТХ = Н.

Это порождает три задачи:

их1 = И, ох 2 = х1, иТх 3 = х2.

Первая система имеет верхнюю треугольную матрицу, вторая — система простых независимых уравнений, третья — нижнюю треугольную матрицу и решается методом прямой подстановки.

Теперь несколько слов нужно сказать об алгоритмах ортогонализа-ции, т. е. алгоритмах, позволяющих построить СЯ-разложение матрицы. Подробно о СЯ-разложениях см. в [4], там же представлены эффективные алгоритмы их построения.

Отражения Хаусхолдера имеют вид:

Р = I - 2 ииТ / иТи,

где Р — произвольная прямоугольная матрица, V — и-мерный вектор.

Удобнее (см. [4]) положить V = х + )||Х|2 е1..

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

Вращения Гивенса имеют вид:

1 ... 0

Є(ї, к ,в) =

0

0

0

0

0 ... -в ... с ... 0

0

і

к

0

0

1

Xi ~Хь

где с = . =, s =

7 2 2 ' I 2 2

Х1 + \Х1 +

Такие преобразования позволяют обнулить выбранный элемент исходной матрицы.

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

В нашем случае ОЯ-разложения используются для:

1) обновлений по времени множителей Холецкого матрицы ковариации неопределенности оценки;

2) обновление по наблюдениям множителей Холецкого информационной матрицы оценки;

3) комбинированных обновлений (по наблюдениям и по времени) множителей Холецкого матрицы ковариации неопределенности оценки.

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

27

Список литературы

1. Мышерин А. В. Применение фильтра Калмана для оценки полного электронного содержания ионосферы по данным GPS наблюдений // Вестник Ка-линингр. гос. ун-та. Сер.: Информатика и телекоммуникации. Калининград: Изд-во КГУ, 2005. Вып. 1-2. С. 140-144.

2. Mohinder S. Grewal, Angus P. Andrews. Kalman Filtering: Theory and Practice. John Wiley & Sons. Inc. 2001.

3. Льюнг Л. Идентификация систем. Теория для пользователя. М.: Наука, 1991.

4. ГолубДж., Ван Лоун Ч. Матричные вычисления. М.: Мир, 1999.

5. Балакришнан. А. В. Теория фильтрации Калмана. М.: Мир, 1988.

6. Eli Brookner. Tracking and Kalman Filtering Made Easy. John Wiley & Sons. Inc. 1998.

7. Mohinder S. Grewal, Lawrence R. Weill, Angus P. Andrews. Global Positioning systems, Internal Navigation, and Inegration. John Wiley & Sons. Inc. 2001.

Об авторах

А. Ф. Лаговский — канд. техн. наук, проф., РГУ им. И. Канта. А. В. Мышерин — асп., РГУ им. И. Канта.

И. И. Шагимуратов — канд. физ.-мат. наук, ЗО ИЗМИРАН.