оптимального оценивания
к задачам навигации
МАТЕМАТИЧЕСКИЕ
ОСНОВЫ
НАВИГАЦИОННЫХ
СИСТЕМ

МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
имени М.В. ЛОМОНОСОВА
А.А. Голован, Н.А. Парусников
МАТЕМАТИЧЕСКИЕ ОСНОВЫ
НАВИГАЦИОННЫХ СИСТЕМ
Часть II
Приложения методов оптимального оценивания
к задачам навигации
2-е издание, исправленное и дополненное
МОСКВА-2012

УДК 527:519.8
ББК 39.471.1:22.18
Г61
Голован А.А., Парусников Н.А.
Г61 Математические основы навигационных систем: Часть II:
Приложения методов оптимального оценивания к задачам на¬
вигации. - 2-е изд., испр. и доп. - М.: МАКС Пресс, 2012. -
172 с.
ISBN 978-5-317-04224-0
Представленные материалы составляют вторую часть книги по мате¬
матическим основам навигационных систем. Она посвящена методам
оптимального оценивания и их приложениям к задачам навигации.
Материалы основаны на учебно-методических разработках кафедры
прикладной механики и управления, лаборатории управления и нави¬
гации механико-математического факультета МГУ имени М.В. Ломо¬
носова в области теории и практики интегрированных навигационных
систем. В этих разработках применен, в свою очередь, опыт многолет¬
него сотрудничества лаборатории и кафедры с рядом ведущих россий¬
ских научно-производственных, научно-исследовательских организа¬
ций, занимающихся проектированием навигационных комплексов.
Для студентов, аспирантов и широкого круга специалистов, зани¬
мающихся прикладными задачами навигации.
УДК 527:519.8
ББК 39.471.1:22.18
ISBN 978-5-317-04224-0
© Голован А.А., Парусников Н.А., 2008
© Голован А.А., Парусников Н.А.,
с изменениями, 2012

Содержание
1	Предисловие	6
2	Задача коррекции в инерциальной	навигации	7
2.1 Общая постановка задачи коррекции . ,	 8
2.2 Задача коррекции и наблюдаемость  	 17
2.3 Критерии наблюдаемости 	 19
3	Методы оптимального и субоптимального оцени¬
вания линейных систем	25
3.1 Элементы теории вероятности и теории случай¬
ных процессов	 26
3.1.1	Понятие случайной величины и ее харак¬
теристики 	 26
3.1.2	Основные законы распределения	 31
3.1.3	Случайные процессы и их характеристики	34
3.1.4	Анализ случайных процессов	 36
3.1.5	Процессы с ортогональными приращения¬
ми. Белый шум	 39
3.2 Стохастические модели линейных динамических
систем	 43
3.2.1	Дискретный случай 	 44
3.2.2	Непрерывный случай 	 44
3.2.3	Дискретизация непрерывных случайных
процессов	 47
3.3 Алгоритмизация задачи оценивания	 50
3.4 Решение переопределенных систем линейных ал¬
гебраических уравнений. Вероятностная интер¬
претация 	 52
3.4.1	Метод наименьших квадратов	 52
3.4.2	Критерий максимального правдоподобия .	53
3.4.3	Задача сглаживания экспериментальных
данных методом МНК при помощи куби¬
ческих сплайнов	 55
3.4.4	Критерий ортогональности и критерий
условного среднего	 58
3.5 Дискретный фильтр Калмана	 62
3.5.1	Алгоритм дискретного фильтра Калмана .	63
3

3.6 Реализация дискретного фильтра Калмана мето¬
дом квадратного корня 	 68
3.6.1 Численная реализация метода квадрат¬
ного корня фильтра Калмана для S-
верхнетреугольной формы	 70
3.7	Непрерывный фильтр Калмана	 74
3.7.1 Представление уравнения Риккати в виде
линейных уравнений большей размерности	77
3.7.2 Некоторые условия устойчивости фильтра
Калмана	 79
3.8	Задача сглаживания (субоптимальный вариант) .	82
3.9	Меры оцениваемости	 87
3.9.1 Сингулярные числа как меры оцениваемости 88
3.9.2 Мера наблюдаемости на основе ортогона-
лизации Грама-Шмидта	 90
3.9.3 Стохастическая мера оцениваемости	....	92
4 Калибровка и выставка инерциальных навигаци¬
онных систем	96
4.1	Калибровка ИНС	 96
4.1.1 Калибровка бескарданной инерциальной
навигационной системы 	 97
4.1.2 Калибровка ИНС с горизонтируемой плат¬
формой 	112
4.2	Задача выставки инерциальных навигационных
систем	120
4.2.1 Приборный трехгранник неподвижно ори¬
ентирован относительно инерциальной си¬
стемы координат 	121
4.2.2 Задача выставки БИНС	125
4.2.3 Выставка ИНС с горизонтируемой плат¬
формой 	142
5 Основные виды навигационной информации, до¬
полнительной по отношению к инерциальной.
Формирование сигналов коррекции	151
5.1	Информация о высоте 	151
5.2	Угловая информация	151
4

5.3	Скоростная информация	154
5.4	Формирование сигналов корректируемых ИНС с
использованием спутниковой навигационной ин¬
формации 	155
5.4.1 Формирование вектора коррекции при по¬
мощи вторичной позиционной спутнико¬
вой информации	 156
5.4.2 Формирование вектора коррекции при по¬
мощи вторичной скоростной спутниковой
информации	 156
5.4.3 Формирование вектора коррекции при по¬
мощи первичной позиционной спутнико¬
вой информации	161
5.4.4 Формирование сигнала коррекции при по¬
мощи первичной скоростной спутниковой
информации	164
5.5	Заключительные замечания	167
5

1.	Предисловие
В пособии излагается теория корректируемых инерциальных
навигационных систем, а также некоторые вопросы калибровки
и выставки инерциальных навигационных систем.
Книга включает описание как методов оптимального оце¬
нивания, традиционно применяемых при решении указанных
выше задач, так и описание конкретных алгоритмов.
Книга является продолжением изданного ранее пособия
по теории инерциальной навигации [1] и основана на учебно¬
методических разработках кафедры прикладной механики и
управления, лаборатории управления и навигации механико¬
математического факультета МГУ имени М.В. Ломоносова.
Авторы выражают благодарность сотруднику лаборатории
Н.Б. Вавиловой за помощь.
С авторами можно связаться по адресу:
E-mail: aagolovan@yandex.ru, www.navlab.ru.
6

2.	Задача коррекции в инерциальной
навигации
Навигационные системы, в основу которых положен метод
инерциальной навигации, обладают рядом достоинств, обеспе¬
чивающих им широкое применение в авиации, судовождении, а
также в других областях, например, при навигации дефектоско¬
пов в газовых и нефтяных трубах, аэрогравиметрии. Наиболее
существенные из них следующие:
• в системах инерциальной навигации формируются прак¬
тически все необходимые параметры для осуществления
управляемого движения объекта — координаты, скорости,
углы ориентации корпуса объекта относительно некоторо¬
го заданного трехгранника;
• работа системы инерциальной навигации не зависит от
внешних условий, сопутствующих движению объекта, та¬
ких, например, как состояние атмосферы, радиопомехи и
т. п.;
• системы инерциальной навигации могут функциониро¬
вать в течение всего времени движения объекта неза¬
висимо от работы других навигационных приборов и
устройств.
С другой стороны инерциальные навигационные системы
(ИНС) имеют ряд недостатков, главные из которых таковы:
• точность определения навигационных параметров в систе¬
мах инерциальной навигации зависит от погрешностей на¬
чальных условий и погрешностей установки элементов си¬
стемы в исходное (начальное) положение;
• ошибки в определении навигационных параметров накап¬
ливаются со временем из-за кинематических ошибок ори¬
ентации приборной системы отсчета, вызванных уходом
гироскопов или инструментальными погрешностями дат¬
чиков угловой скорости (ДУС).
7

Именно эти недостатки привели к тому, что в системах инер¬
циальной навигации для улучшения их динамических и точ¬
ностных свойств стали использовать помимо инерциальной ин¬
формации дополнительную информацию неинерциальной при¬
роды. Такой дополнительной информацией являются: позици¬
онная, доставляемая радиосистемами различного типа; ско¬
ростная, доставляемая доплеровскими измерителями скорости
(ДИС), лагами, одометрами; угловая, которую доставляют, на¬
пример, астросистемы; информация о высоте.
В настоящее время получили самое широкое распростра¬
нение спутниковые навигационные системы (СНС), доставля¬
ющие позиционную и скоростную информацию. Существуют
модификации приемников сигналов СНС, позволяющие с по¬
мощью разнесенных антенн определить угловую ориентацию
объекта относительно инерциального пространства. СНС стали
основным источником дополнительной информации для кор¬
рекции ИНС. Спутниковым навигационным системам посвящен
ниже специальный раздел.
2.1.	Общая постановка задачи коррекции
Рассмотрим постановку задачи коррекции ИНС в общем виде
[2]. Пусть поведение динамического объекта подчиняется урав¬
нению
X = F(X,U),	(2.1)
где X - вектор состояния, U - внешнее воздействие, доступное
измерению.
Для определения текущего значения вектора состояния X
привлекаются два вида информации, условно называемые ос¬
новной и дополнительной:
1.	Основная информация
X' = X{t0) + a*, U'{t) = U(t) + u(t),
где хо - ошибка начального знания вектора состояния, u(t)
- инструментальная погрешность измерения величины U.

2. Дополнительная информация
Z*(t) = e(X(t)) + r(t),
где r(t) - инструментальная погрешность.
Принципиально важно, что совокупная информация Xq, U',
Z* избыточна, где избыточность понимается как неоднознач¬
ность определения вектора состояния при идеальной исходной
информации
хо = 0, u(t) = 0, r(t) = 0.
Вопрос состоит в том, как построить систему, определяю¬
щую вектор состояния с наилучшей в каком-то смысле точно¬
стью. В [1] в связи с задачей привлечения информации о высоте
при построении алгоритмов ИНС рассмотрены три схемы обра¬
ботки совокупной информации. В первой - за счет дополни¬
тельной информации меняется функциональная структура си¬
стемы, во второй - часть основной информации заменяется до¬
полнительной. За счет последней в обоих случаях меняются в
лучшую сторону динамические свойства системы, но избыточ¬
ность информации полностью не используется.
Рассмотрена также третья схема, когда в модельные уравне¬
ния ИНС вводится обратная связь, формируемая при помощи
дополнительной информации. За счет этой связи улучшаются,
как и в первых двух схемах, динамические свойства системы.
Но в рамках третьей схемы можно поставить задачу построе¬
ния оптимального алгоритма, позволяющего минимизировать
ошибки определения навигационных параметров. Это означает,
что возможности, заложенные в совокупной информации, реа¬
лизуются полностью. При этом погрешности хо, u(t), r(t) пред¬
полагаются настолько малыми, что задача оптимизации может
быть решена в линейной постановке.
Структуры систем, информационно эквивалентных друг
другу, могут быть различны. Одна из них состоит в следую¬
щем. Пусть основная информация вводится в систему, которую
будем обозначать /, в соответствии с модельными уравнениями
X' = F(X',U'),	X'(t0)	=	X'0.	(2.2)
9

Уравнения ошибок такой системы относительно вектора х(£) =
X'(t) — X(t) имеют вид
х = Ах + q	(2.3)
где
,	dF(X',U')	dF(X',Uf)
А=  gxi .	4=	off,	u> 4to)=x0.
Введем величину z = Z* — Q(X'). Разлагая величину
в(Х') = 0(Х + х) в ряд Тейлора и оставляя только линейные
члены разложения, получим
z = Нх + г,	(2.4)
где
двщ^дерп
эх дхг '
Задача сводится теперь к определению (оцениванию) вели¬
чины х, удовлетворяющей уравнению (2.3) при помощи вели¬
чины z (2.4).
Пусть выбраны математические модели инструментальных
погрешностей q и г, соответствующий этим моделям крите¬
рий, и построен алгоритм (оператор) L, доставляющий оценку
х = L[z\ такой, что он минимизирует ошибку оценки Ах = х —х
с точки зрения выбранного критерия. Тогда оценка X вектора
состояния находится из соотношения
X = X' - х = X + Ах.	(2.5)
Ошибка Ах оценки величины х оказывается, таким образом,
ошибкой оценки вектора состояния X.
Величину z будем называть вектором коррекции, а саму за¬
дачу определения величин х, X - задачей коррекции.
Замечание. Представляется очевидным, что любое реше¬
ние по определению вектора состояния X при помощи сово¬
купной информации Xq, Uf, Z* в рамках принятых моделей не
улучшаемо по сравнению с только что описанной процедурой.
10

Тем не менее, у некоторых специалистов время от времени
возникает идея улучшить точность системы за счет введе¬
ния в систему I обратных связей, образованных с использова¬
нием дополнительной информации Z*.
В триаде "математическая модель инструментальных
погрешностей", "критерий", "алгоритм" первичной является
модель погрешностей. Пока не сформирована такая модель, за¬
дача определения вектора состояния не может быть коррект¬
но поставлена. Описание такой модели в инженерной практике
называют обычно техническими условиями. Выбор адекватной
модели погрешностей составляет наиболее трудную и нефор¬
мальную часть задачи.
Обычно конкретный физический смысл составляющих эле¬
ментов модели погрешностей понятен разработчику прибора,
но это часто не принимается во внимание, а модель погрешно¬
стей строится из здравого смысла и предшествующего опыта,
то есть носит феноменологический характер.
Как правило, в качестве математической модели инструмен¬
тальных погрешностей выбираются линейные дифференциаль¬
ные уравнения с известными коэффициентами. Те составляю¬
щие погрешностей, которые не укладываются в такую модель,
входят в уравнения как возмущения. Эти возмущения будем
называть немоделируемыми.
Итак, принимается такая модель величин q и г, что уравне¬
ния (2.3) и (2.4) приводятся к виду:
х = Ах + Tiaei + £1,
aei = r2aei+C2,
эе2 = Г3ае2 + Сз,
z = Нх Н- Г4зе2 -Ь £.	(2*6)
Здесь Г» - известные матрицы, в общем случае зависящие от
времени; aei, ае2 — параметры моделируемых составляющих ве¬
личин q, г; составляющие £г (г = 1,2,3), е — немоделируемые
возмущения, относительно которых должны быть приняты те
или иные гипотезы.
Чаще всего полагается, что эти величины <д, £2? £3, е - белые
шумы с известными интенсивностями. О понятии белого шума
11

будет рассказано в следующем разделе.
Введенная соотношениями (2.6) модель инструментальных
погрешностей достаточно универсальна и гибка и, как показы¬
вает опыт проектирования навигационных систем, весьма кон¬
структивна. Ее частным случаем является широко использу¬
емая модель, в силу которой компоненты инструментальных
погрешностей представляются суммой двух слагаемых: первое
слагаемое - конечномерное разложение по известным базовым
функциям, второе слагаемое - стационарный случайный про¬
цесс с дробно-рациональной спектральной плотностью. Описан¬
ные в [1] математические модели инструментальных погрешно¬
стей ИНС имеют именно такую структуру.
Обобщением указанной модели служит модель, в кото¬
рой некоторые элементы матриц Г* неизвестны. Такая модель
при построении алгоритмов навигационных систем приводит к
адаптивным задачам.
Иная модель предполагает введение ограничений на компо¬
ненты векторов q и г (а также, может быть, на их первые , вто¬
рые и т.д. производные). Подобная модель обычно связывается
с минимаксным критерием и задачей о накоплении возмуще¬
ний. Две последние модели играют при построении алгоритмов
оценки чаще всего вспомогательную роль.
Введем обозначения
Задача состоит теперь в том, чтобы на основании (2.8) постро¬
ить наилучшую в некотором смысле оценку £ вектора £ и тем
самым получить оценку х для вектора х.
В качестве алгоритма L, доставляющего оценку £ вектору
£, в большинстве случаев выбирается оператор, описываемый
следующими соотношениями
Щ = (Я ог4).
В этих обозначениях перепишем соотношения (2.6)
(2.7)
£ = А^ + С, z = Htf + e.
(2.8)
^ = A^ + K(z-H^),	*(*>)	=	6>.	(2-9)
12

Уравнение ошибок оценки относительно величины Д£ = £ — £
имеет вид
Д£ = (А* - КЩ)А£ -Ке + С	(2.10)
Матрица К называется матрицей усиления или коэффици¬
ентом усиления. Она подлежит выбору из условий минимиза¬
ции, в том или ином смысле, величины Д£. О достоинствах это¬
го алгоритма речь пойдет ниже. Пока же отметим его очевидное
свойство несмещенности — тождественное равенство
Д£ = 0,
при
С = 0, 8 = 0, £(*<,) = £(to).
Итак, сформулирована задача коррекции, возникающая при
оценивании состояния динамического объекта, когда привлека¬
емая для этой цели информация избыточна. Описан вариант ре¬
шения задачи, состоящий в том, что с помощью дополнительной
информации оценивается выходная ошибка системы J, входом
которой служит основная информация.
Однако возможны иные формы приборного решения ука¬
занной задачи, информационно полностью эквивалентные рас¬
смотренному варианту - формы, в которых используется та же
исходная информация, и вектор состояния динамического объ¬
екта определяется с той же точностью.
Речь идет об использовании, хотя бы частичном, интеграто¬
ров основной системы I для интегрирования дифференциаль¬
ных уравнений оценок. Это означает, что в систему вводятся
дополнительные (обратные) связи и дополнительная информа¬
ция является также входом системы.
Иначе говоря, в подобных системах задача коррекции ре¬
шается не просто как задача оценивания, но и как задача, при
решении которой используются обратные связи, вводимые в си¬
стему [2], [3]. Эти связи зависят от дополнительной информа¬
ции. Прежде чем описывать такие системы, запишем более по¬
дробно соотношения (2.9), (2.10).
Уравнения коррекции.
х = Ах + Tiaei + K\(z — Нх — Т^2)^
13

aei = r2aei + K2{z - Hx - Г4аё2),
ae2 — Г3ае2 + Кз{г — Hx — Г4аэ2).	(2*11)
Уравнения ошибок оценки.
Ад; = ААх + TiAaei — Ki(HAх + Г4Дае2) — К\г + £i,
Aaei	= r2Aaei - K2(HAx + Г4Аае2) - K2e + C2,
Дае2	= ГзДае2 — К$(НАх + Г4Дае2) — К$е + Сз-	(2.12)
Рассмотрим теперь систему /*, в которой в блок I введены
обратные связи. Пусть алгоритм системы I* описывается мо¬
дельными уравнениями
X* = F(X\U') - Ki(z* - Г4ае2) - Тг&и
аэ 1	=	r2aei 4- K2(z* — Г4ае2),
зе2 = Г3аё2 + K%(z* — Г4ж2),
z* = Z* — в{Х*) = Нх* + Г4ае2 + е,	(2.13)
где х* =	X — X* - ошибка определения системой I*	вектора
состояния X.
Уравнения ошибок системы I* относительно величин
ж*, Aaei, Дзе2 имеют вид
х*	=	Ах* + TiAaei - Кг{Нх* + Г4Дае2) - Кхе	+	Сь
Aaei	=	r2Aaei - К2(Нх* + Г4Дае2) - К2е + £2,
Дае2	=	Г3Дае2 — К2(Нх* + Г4Дае2) — К^е + Сз-	(2-14)
Уравнения (2.12), (2.14) структурно полностью совпадают при
должном задании начальных условий
х* — Ах.
Таким образом, системы Г и I* информационно эквивалент¬
ны, то есть обеспечивают одинаковую точность решения задачи.
При этом число дифференциальных уравнений, описывающих
работу системы /*, меньше числа уравнений, соответствующих
системе на величину размерности вектора X.
Возможны иные способы построения алгоритма системы, в
которой используется избыточная информация, например, при
14

введении обратных связей по части составляющих вектора X.
Эквивалентные обратные связи могут вводится также в случае,
когда в уравнении (2.1) используются избыточные переменные,
между которыми существуют алгебраические связи. Примером
такого вида уравнений служат, например, матричные кинема¬
тические уравнения Пуассона.
Выделим один из вариантов, достаточно часто встречаю¬
щихся в навигационной практике, когда при помощи вектора
коррекции в систему I вводятся обратные связи с целью сде¬
лать уравнения ошибок асимптотически устойчивыми и тем са¬
мым оправдать гипотезу линейности. Далее задача минимиза¬
ции ошибок оценки решается по уже описанному выше образцу.
Пусть модельные уравнения системы имеют вид
X* = F(X*, U') + Cz\ z* = Z - 0(Х*) = Нх* + г. (2.15)
Уравнения ошибок относительно величины х* = X* — X имеют
вид
х* = (А + СН)х* + q + Сг.	(2.16)
Здесь величина С выбирается из условия улучшения динамиче¬
ских свойств уравнений ошибок. Далее, как и при выборе опера¬
тора L, может быть поставлена задача получения оптимальной
оценки х* вектора я* при помощи вектора коррекции z*.
Поскольку задача коррекции может быть решена структур¬
но различным образом, встает вопрос о выборе подходящей
структуры. Но такой выбор не может быть сделан, исходя из
формальных линейных критериев. Здесь могут играть роль
традиции, соображения надежности, взаимоотношения коллек¬
тивов, проектирующих отдельные блоки системы и т.д. В по¬
следние годы, когда возможности бортовых вычислителей зна¬
чительно выросли, для инерциальных навигационных систем
с прецизионными чувствительными элементами задачу коррек¬
ции чаще всего предпочитают решать как задачу чистого оцени¬
вания. Но возможны варианты, когда ошибка xq не достаточно
мала или не достаточно малы инструментальные погрешности
при использовании грубых инерциальных датчиков, таких как
MEMS (MicroElectroMechanical Systems) датчиков. Тогда ока¬
15

зывается предпочтительнее варианты алгоритмов с введением
обратных связей.
Замечания:
1.	Если сигнал коррекции z не скаляр, тпо вообще гово¬
ря, существует теоретическая возможность декомпо¬
зиции (расщепления) алгоритма коррекции на ряд более
простых алгоритмов, число которых равно размерности
корректирующего вектора (см., например, Ц]). Но такую
декомпозицию вряд ли стоит использовать за исключе¬
нием случая, когда она явно напрашивается. Пример та¬
кой декомпозиции будет описан ниже при обсуждении
задачи начальной выставки ИНС. Опишем примененный
там прием в общем виде. Пусть структура уравнений,
описывающих задачу коррекции имеет вид
х/	=	A\Xi + CihfjXu + qi,
xi i	=	А2хц + C2hjxi + gj/,
Z\	=	/i/Xj+ri,
z2	=	hfjXn + r2.
Субоптимальный алгоритм может быть построен в
следующем виде
xi = Axxi + C\z2 + Ki(zi - hjxi),
хц = А2хц + C2z\ 4- K2(z2 — hJjXn).
Соответствующие уравнения ошибок оценки имеют вид:
Дх/ = (Ai - Kihf )Дх/ + C\r2 + qi - K\r^
Дх/j = (А2 — K2hjj)Axn + C2r\ + qn — K2r2.
Здесь К\ и К2 выбираются из условия оптимальности
алгоритмов оценивания двух подсистем, которые пола¬
гаются независимыми друг от друга. Эти системы, оче¬
видно, таковы:
X/ = Aix} + qi + C\r2, zl = hjx} + ri,
X//	=	A2x^j	+ qn + C2r\,	= hfjx’jir2.
16

2. Несколько слов о задании матриц интенсивности белых
шумов, входящих в уравнения ошибок ИНС и вектор кор¬
рекции. Априори они бывают в какой-то степени извест¬
ны. Но лучше подобрать их в результате моделирования
задачи коррекции с использованием реальных измерений.
2.2.	Задача коррекции и наблюдаемость
При решении коррекционной задачи, особенно при интерпрета¬
ции результатов работы соответствующих алгоритмов, оказы¬
вается полезным понятие наблюдаемости. Прежде чем форма¬
лизовать это понятие, поясним его на простом примере.
Пусть требуется построить оценку вектора х = (х\,Х2)Т,
удовлетворяющего уравнению:
х\ = Ai#i , Х2 = А2X2
при помощи измерения
z(t) = xi (t)+x2(t).
Для определения х достаточно найти начальные условия
яиь#20- Тогда
xi(t) = xi0eXlt ,	x2(t)	=	X2oeX2t.
Для вычисления яю,#20 имеем соотношение
z(t) = xi0eXlt + x2oeX2t.
Самый простой способ вычисления этих величин использовать
измерение z(t) в два различных момента времени, например,
при t = 0 и t = 1. Тогда имеем
z(0) = хю + х2о , z( 1) = xi0eAl + х2оеЛ2.
Понятно, что хю и Х20 определяются однозначно, если Ai ^ А2.
Если же Ai = А2, то можно определить только сумму хю+я2о, а,
стало быть, только сумму Xi(t)+X2{t). В первом случае система
наблюдаема, во втором - ненаблюдаема.
17

Если величины Ai, А2 близки:
A2 = Ai(l + e), (е — малый параметр),
то хотя формально имеет место наблюдаемость, задача плохо
обусловлена и можно говорить о слабой (плохой) наблюдаемо¬
сти фазовых переменных задачи.
Далее понадобятся некоторые сведения из теории линейных
дифференциальных уравнений.
Пусть вектор х удовлетворяет уравнению
х = A(t)x + (p(t),	(2-17)
с начальными условиями x(to) = #о*
Решение уравнения (2.17) имеет вид
x(t) = x^t) + xn(t),
где xj(t) - решение однородного уравнения
xj = A(t)xi,
с заданными начальными условиями #/(£о) = #о-
Оно имеет вид
xj(t) = Ф(^о)хо-
Матрица Ф(£,£о) носит название переходной матрицы и удовле¬
творяет уравнению
дФ^1о) = A№(t, to), Ф(*0, to) = Е.
Некоторые свойства переходной матрицы:
1. det Ф(г,£0) ф 0;
2. Ф(£,т)Ф(т, s) = Ф(£, 5), откуда следует Ф(s,t) = Ф-1(^, s);
3 дФ£й=_Атт{зЛ
18

При А = const переходная матрица Ф(£До) имеет вид
Ф(Мо) = ел<*"4
где eAt - матричная экспонента, определяемая рядом
Вектор xu(t) удовлетворяет неоднородному уравнению
хп = A(t) + <p(t),
с нулевыми начальными условиями. Его решение имеет вид
xn{t) = Ф (t,r)(f{r)dT.
2.3.	Критерии наблюдаемости
Далее будет сформулировано понятие наблюдаемости, сформу¬
лированы и доказаны критерии наблюдаемости и указаны неко¬
торые свойства наблюдаемых систем [4]. Доказательства здесь
приводятся потому, что они полезны для понимания теории кал-
мановской фильтрации и сглаживания.
Итак, рассматривается система, поведение которой подчи¬
няется уравнению
Система (2.18) называется вполне наблюдаемой в момент t, если
существует конечный момент to < t такой, что можно однознач¬
но определить состояние системы x(t) из наблюдения выходной
функции z(t) на отрезке [to,t\.
В определении присутствует слово "вполне1’. Оно введено
для того, чтобы отличить качество системы от того факта, что
в систему введено наблюдение z(t), то есть в системе произво¬
дится наблюдение. При общении специалистов такие тонкости
не учитываются и слово "вполне” опускается, и это никогда не
приводит к недоразумению.
x(t) = A(t)x(t) , z(t) = H(t)x{t).
(2.18)
19

Поскольку свойство наблюдаемости целиком определяется
свойством пары (A, if), то вместо того, чтобы сказать, что си¬
стема (2.18) вполне наблюдаема, говорят, что наблюдаема пара
(ЛЯ).
Сформулируем критерий наблюдаемости. Образуем сим¬
метрическую матрицу Af(t,to) > 0, которая называется грами-
аном наблюдаемости:
t
•W(Mo) = Уфт(т,£)#т(т)#(т)Ф(т,*)йт.	(2.19)
to
Из (2.19) следует уравнение, которому удовлетворяет грами-
ан наблюдаемости
9M^to) = -AT(t)JV(t,to)-W,t0)A(t)+HT(t)H(t), Лf(t,t) = 0.
(2.20)
Имеет место свойство инвариантности ранга матрицы N при
невырожденном стационарном преобразовании вектора состоя-'
ния.
Положим £ = Сх, С = const, det С Ф 0. Очевидны уравнения
£ = С АС_1£, г = ЯС"1^	(2.21)
Грамиан системы (2.21) обозначим через ЛГ*. Используя
(2.19) и (2.21), получим
ф* = сфст, к = (с_1)тл/’с'_1, ф*л/; = сфлгс-1. (2.22)
Таким образом, N и ЛГ* конгруэнтны. Но конгруэнтные матри¬
цы имеют одинаковый ранг.
Теорема. Для наблюдаемости пары (А, Я) в момент t необ¬
ходимо и достаточно, чтобы нашелся момент to < t такой, что
detAf(t, t0)	0,	т.е.	Af(t,t0) > 0.
Доказательство. Рассмотрим уравнения (2.18)
х = Ах, z = Ях.
Зафиксируем некоторый момент to как начальный. Для опре¬
деления x(t) имеем уравнение
z(T) = H(r)$(r,t)i(t),	(2.23)
20

которое должно удовлетворяться для всякого т из интервала
Достаточность. Попробуем построить конкретный алгоритм
определения x(t) из уравнения на интервале [*оД]- Выберем, в
качестве пробной, некоторую величину x(t). Для произвольно¬
го момента т 6 [to,t] это пробное значение x(t) не обязатель¬
но удовлетворяет уравнению (2.23) и имеется рассогласование
(невязка)
Естественно потребовать, чтобы пробное решение x(t) миними¬
зировало в некотором смысле совокупность всех невязок на ин¬
тервале [£о> t]. Мерой невязки может служить следующий функ¬
ционал
Полученное нами решение (2.23) является точным. Это следует
из прямой подстановки выражения (2.26) в выражение (2.23).
Отсюда также следует, что наш функционал на решении (2.26)
обращается в ноль: J = 0.
Необходимость. Предположим, что система наблюдаема, т.е.
существует единственное решение уравнения (2.23)
Пусть матрица ЛДМо) не является положительно определен¬
ной, т.е. для некоторого х
[to,*].
Az(t) = z(t) — #(т)Ф(т, £)ж(£).
t	t
J[x] = J \\Az(r)\\2 dr = J Azt(t)Az(t) dr. (2.24)
to
Решая задачу безусловной минимизации, имеем
(2.25)
которое приводит при det М ^ 0 к решению
t
x(t)=J\f 1{t,t0) J Фт(т, t)HT (t)z(t) dr. (2.26)
to
z(t) = Я(т)Ф(т,*0)я(*), т € [to,*].
xTAf(tjto)x = 0,
(2.27)
21

и этому х отвечает свое измерение z(r) в соответствии с соот¬
ношением (2.23). Из (2.27) следует:
t	t
J	{r,t)HT{т)Н(т)Ф{т^)дт	х = J zT(t)z(t) dr = О,
_io	-	*o
или z(r) = 0 для всех т Е [to,t\. Следовательно, для рассматри¬
ваемой системы существует решение x(t), которое невозможно
отличить от нулевого, что противоречит условию наблюдаемо¬
сти.
Замечание. Вместо уравнения (2.23), определяющего x(t),
можно рассматривать уравнение для определения x(to) = Xq
z(t) = Я(т)Ф(т,to)x(to), to <Т <t.
Эквивалентный грамиан наблюдаемости выглядит так
t
■A/’tMo) = J ^T(r, to)HT (т)Н(т)Ф(т, to) dr.	(2.28)
Теорема.	Для наблюдаемости пары (А, Н)	в момент t необ¬
ходимо и достаточно, чтобы нашелся момент to < t такой, для
которого уравнение относительно ае:
Я(т)Ф(т, t)se = 0, to<r<t,	(2.29)
имело бы при всех т единственное тривиальное решение.
Поясним информационный смысл грамиана наблюдаемости.
Он состоит в том, что грамиан позволяет судить не только о
принципиальной возможности осуществлять наблюдение за со¬
стоянием системы, но и, до некоторой степени, о качестве на¬
блюдения. Рассмотрим величину
t
J = J zT(r)z(r)dr,	(2.30)
которую можно интерпретировать как "энергию” выходного
сигнала, накопленную на интервале [to,t]. Имеет место очевид¬
ное соотношение J = xT(t)Af(t, to)x(t).
22

Существует ортогональное преобразование ае = Fx(t) такое,
что
где A j — корни характеристического полинома det(A Е — М).
При малом значении A j вклад координаты ае^* в энергию вы¬
ходного наблюдаемого сигнала также мал, и следует ожидать,
что ошибка в определении этой координаты, порожденная по¬
грешностью измерения z, будет обратно пропорциональна соб¬
ственному числу A j, т.е. при малом A j координата ае^ - слабо
наблюдаемая величина.
Критерий наблюдаемости в стационарном случае.
Теорема. Для наблюдаемости стационарной пары (А, Н) необ¬
ходимо и достаточно выполнения условия:
Матрица N называется матрицей наблюдаемости.
Критерий поясним рассуждениями, которые, по существу,
представляют вариант его вывода. Аналитическая функция z(t)
однозначно определяется значениями этой функции и всех ее
производных в некоторый момент t. В силу уравнения (2.18),
получим в этот момент
Из теории матриц известна теорема Гамильтона-Кэли, которая
гласит:
любая квадратная матрица удовлетворяет своему собственному
характеристическому уравнению.
Это означает следующее. Запишем характеристическое
уравнение матрицы А
п
н \
rank N = n, N =
НА
\ НАп~1
z = Нх, i = ЯАх,..., z^ = НА^х, ....	(2.31)
det(XE — А) — А71 -Ь ai\n * + ... -Ь ап = 0.
23

Имеет место соотношение
Ап -I- aiAn ^	+	апЕ	=	0.
Тогда все производные z начиная с	в	силу	теоремы
Гамильтона-Кэли оказываются линейными комбинациями ве¬
личин z, z7, zn, ...z^n~^\ Поэтому бесконечная система уравне¬
ний (2.31) относительно х равносильна конечной системе
z = Hx, z = НАх,... z(n_1) = НА^х. (2.32)
Таким образом, возможность определения вектора х по изме¬
рениям z(t),t £ [£о,оо) эквивалентна возможностям, содержа¬
щимся в системе (2.32). Но необходимое и достаточное условие
для того, чтобы система (2.32) имела единственное решение та¬
ково:
rank N = dim(x).
Основное свойство стационарной системы с наблю¬
даемой парой (А, Н). Оценку х вектора х будем искать при
помощи динамического несмещенного алгоритма оценивания
х = Ах + К (z — Нх),	(2.33)
о котором уже говорилось ранее.
Уравнение ошибок оценки относительно величины Ах = х —
х имеет вид
Ах —{А — КН)Ах.
Если пара (А, Н) - наблюдаема, то коэффициенты характери¬
стического уравнения
det(XE-(A-KH))=0
за счет выбора элементов матрицы К могут быть сделаны лю¬
быми наперед заданными. В частности, их можно выбрать так,
чтобы корни указанного характеристического уравнения име¬
ли отрицательные действительные части. Но тогда уравнение
ошибок оценки асимптотически устойчиво, то есть Ax(t) —» 0
при t —> оо. Доказательство этого факта содержится в [4].
24

3.	Методы оптимального
и субоптимального оценивания
линейных систем
Раздел содержит описания методов оценивания, широко ис¬
пользуемых при решении многих прикладных задач навигации,
таких как коррекция ИНС при помощи информации неинер-
циальной природы, выставка ИНС, калибровка инструменталь¬
ных погрешностей, спутниковая навигация.
Поскольку эти методы основаны на гипотезе о стохастиче¬
ской природе инструментальных погрешностей, в разделе при¬
водятся минимально необходимые сведения из теории вероят¬
ности и теории случайных процессов. В конечном счете речь
здесь идет об оценивании вектора состояния ж, связанного либо
переопределенной системой алгебраических уравнений
z = Нх + г,	(3.34)
либо, в более сложном случае, динамическими уравнениями
х = Ах + <7,	(3.35)
и измерениями вида (3.34).
Дискретный вариант:
1) = *^jx{tj) Qji z(tj) = Hjx(tj)	-f- Tj.	(3.36)
Излагаемым	ниже	материалам посвящены	многие пособия.
Среди них следует выделить работу [5], к сожалению, не пере¬
веденную на русский язык.
Здесь мы следуем работе [4]. Предлагаемый текст - резуль¬
тат переработки двух глав из указанной книги. В раздел не
вошли частотные методы построения алгоритмов оценивания и
их анализа. В частности, не описаны такие важные для прило¬
жений вещи, как явления маскировки частот, явления Гиббса,
частотные способы построения формирующих фильтров и т.д.
Возможно, по этим вопросам будет написано отдельное посо¬
бие.
25

3.1.	Элементы теории вероятности и теории
случайных процессов
При строгом изложении теории вероятности вводится понятие
вероятностного пространства вместе с сопутствующей этому
пространству аксиоматикой. Здесь мы следуем упрощенному
пути и первичным считается понятие случайной величины.
3.1.1.	Понятие случайной величины
и ее характеристики
Универсальной характеристикой, описывающей случайную ве¬
личину X, служит функция распределения F(x). По определе¬
нию,
F(x) = P(X <х),
где Р — обозначение вероятности.
В более подробной записи, когда X = (Xi, Х2,..., Хп)т,
F(x)	F{x\, Х2, •.., хп^	Р{Х-\ <С Х\, Х2	Х2) • • •) Хп <с Хп}.
Случайная величина называется непрерывной, если непре¬
рывна ее функция распределения. Очевидно, что F(x) монотон¬
но не убывающая функция аргумента х.
Если существует функция f(x) = f(x 1, £2,..., жп), такая,
что
Х\ Х2	хп
F{x) = //-/ f (v>i, и2,... ,un)du\du2 .. •dun,	(3.37)
—00 —00	—00
при любом х, то эта функция называется плотностью вероят¬
ности.
Замечание. Здесь и далее записи F и f служат не для
обозначения каких либо конкретных функций, а как символы
функций распределения и плотности.
Из (3.37) следует связь функций распределения и плотности
вероятности
//дЛ   dF (Х\ , Х2 ) • • • 5 %ть)
dx\dx2.-.dxn
26

Переход от многомерных характеристик F(x) и f(x) к одно¬
мерным F(xj), f(xj) (Xj — соответствующая компонента векто¬
ра х = (xi,X2,.. • ,жп)) осуществляется по следующим форму¬
лам (для простоты положено п = 2, j = 1):
оо
F(xi) = F(xi,+oo), f(xi) = J f(xi,X2)dx2.
— OO
Пусть X и У — две случайные векторные величины. Функ¬
ция
F(x/Y = у) = Р(Х < x/Y = у),
называется условной функцией распределения вероятностей X
при условии У.
Определение условной плотности вероятности.
Имеет место формула Байеса
f(x/v) = ^х’У)
п ,у) № ’
аналогичная формуле, определяющей условную вероятность со¬
бытий.
Если f(x,y) = f(x)f(y), случайные величины X и У на¬
зываются независимыми. Очевидно, для независимых величин
справедливы соотношения
f(x/y) = f(x), f{y/x) = f(y).
Основные числовые характеристики случайной (в общем
случае векторной) величины — математическое ожидание и ко¬
вариация. Пусть д(Х) — некоторая функция случайного (в об¬
щем случае векторного) аргумента X. Математическим ожи¬
данием М[д(Х)\ = д,д функции д(Х) называется следующая
величина:
М\д{Х)\= J 9(x)dF(x),
(ае)
где (ае) — область изменения случайной величины.
27

Если определена плотность вероятности /(х), то
оо оо	оо
оо
—оо —оо	—оо
Величина М [.X] носит название математического ожидания
случайной величины X (другой термин — среднее значение X).
Для того чтобы выделить операцию определения математи¬
ческого ожидания, будем использовать запись М [X], результат
операции обозначим через //х.
Оператор М[-\ линеен, то есть имеет место соотношение
Введем обозначение X = X — //х. Величину X будем назы¬
вать центрированной случайной величиной.
Матрицей ковариации случайного вектора X (или просто
ковариацией) называется матрица Rx:
Эта матрица по определению симметричная и неотрицательно
определенная. В том случае, когда Rx — положительно опреде¬
ленная матрица, вектор X будем называть невырожденным (и
вырожденным в противном случае).
сиями соответствующих компонент и обозначаются Di = D [Xi].
Величины &i = \fDi называются среднеквадратичными
(или стандартными) отклонениями.
Из определения дисперсии следует, что она служит мерой
отклонения скалярной случайной величины относительно свое¬
го среднего значения.
Чтобы выяснить смысл недиагональных элементов матри¬
цы Дх, которые называются моментами корреляции между
М [ад(Х) + bh(Y)] = аМ [д{Х)\ + ЪМ [Л(У)].
В частности,
М [аХ + bY] = аМ [X] + ЪМ [У].
о
о
д* = м[х-хт].
Диагональные элементы
называются диспер-
28

компонентами вектора, рассмотрим двумерный вектор X =
о о
(Хх,Х2)т. Величина R12 = М[Х\Х2] — момент корреляции
между случайными величинами Х\ и Х2.
Одновременно с моментом корреляции R12 используется
безразмерная характеристика г 12, называемая коэффициентом
корреляции. Она получается из момента корреляции нормиров¬
кой:
г 12 =
Д12
<7l<72 ’
2 2
где (^ = М[Х1],	о22=М[Х2].
Заметим, что матрица ковариации Rx может быть представ¬
лена в виде
Rx = & X 1
где
(
<71
0 .
.. 0
\
/ 1
гг 2 •
г1п
\
Ьх —
0
<т2 •
.. 0
II
н
Г\2
1 .
Г2п
1
0
0 .
(Tn
)
\ ггп
Г2п
.. 1
/
Легко видеть, что R\2 (соответственно г 12) = 0, если Х\ и
Х2 — независимые случайные величины.
Замечание. Из некоррелированности не следует независи¬
мость случайных величин. Рассмотрим пример. Пусть ком¬
поненты случайного вектора Y = (Yi,Y2)T являются коорди¬
натами точки Y на плоскости OY1Y2. Будем считать, что
точка Y всегда принадлежит кругу единичного радиуса, центр
которого совпадает с началом координат. Предположим, что
двумерная случайная величина Y = (Yi, Y2)T распределена рав¬
номерно в указанной области. Тогда
Му = 0,	/(2/1,2/2)
7Г
R12 = М \Y iY2\ = J —yiy2dy\dy2 = О,
y?+v5<i
29

поскольку подынтегральное выражение является нечетной
функцией своих аргументов, а область интегрирования цен¬
трально симметрична. С другой стороны, величины yi,y2,
очевидно, зависимы, поскольку
У\+У2< !•
Попытаемся Хч выразить через Х\ линейным образом:
Хч ~ аХ 1 + 6,
подобрав коэффициенты а и 6 так, чтобы среднеквадратичная
ошибка S линейного представления была минимальной:
82 — М [(Хч — аХ 1 — 6)2] = min.
ayb
Получим
S2 = <т|(1 - г22) + (acr 1 - Г12<т2)2 + (у.2 ~ ащ - Ь)2.
Минимальное значение 82 обеспечивается при
ас 1 — 7*1202 = 0,	Ц2 —	-6 = 0.
И тогда
52 = *2(1-г22).
Ошибка 82 максимальна при 7*12 = 0 (R12 = 0) и равна ну¬
лю при |т*121 = 1. Кроме того, |ri21 < 1, поскольку 82 > 0, и,
соответственно,
|i?l2| < (Т1СТЧ-
Обратим внимание на то, что условие рч — api —6 = 0 экви¬
валентно условию М [Хч — аХ\ — 6] = 0, то есть величина 82
служит дисперсией D от Д, где Д = Х4 — аХi — 6. Условие
М [Д] = 0 называется условием несмещенности. Об этом более
подробно будет сказано далее.
С другой стороны, если имеет место точное равенство
Хч = аХ 1 + 6,
30

то, как легко показать,
т 12 = sign а.
Таким образом, коэффициент корреляции служит мерой ли¬
нейной связи Х2 с Х\. Так как г 12 отражает степень только ли¬
нейной связи, из некоррелированности Х2 и Xi еще не следует
их независимость.
Для комплекснозначного случайного вектора X матрица ко¬
вариации Rx определяется соотношением
Здесь символ * обозначает две операции: транспонирование и
сопряжение.
3.1.2.	Основные законы распределения
Закон равной плотности (для скалярной случайной величи¬
ны). Закон определяется двумя параметрами х\ и Х2 (Х2 > х\):
Нормальный закон распределения (закон Гаусса).
Одномерный случай.
Закон определяется параметрами рх = М [X] и = D [X]:
RX = M XX .
/(*)
/(*)
1
если г 6 [11,12])
если x^[xi,X2],
Ж2 — Xi
О,
Х2 +Х1
2

Для функции Ф(£) составлены подробные таблицы.
В приложениях часто используются два значения вероятно¬
сти, связанные с нормальным законом распределения:
1. Р{\Х-рх\ <ах} = 0.67,
2. Р {\Х — fix\ < Зсгх} = 0.997.
Условие 2 инженерами называется правилом 3<т. Согласно это¬
му правилу пренебрегают вероятностью того, что
\Х - рх| > 3ах,
то есть значения х = рх ± 3ах считаются предельно допустимы¬
ми.
Многомерный случай.
Х = (ХъХ2,...,Хп)т,
f (^) f (*^1 j %2) • • • j *^п)
(27r)_5(det Дх)"5 -ехр| - ^(х - fix)TR^1 (х - цх) ,
цх=М[Х], Rx = M[X■ Хт].	(3.39)
Замечание. Важность нормального распределения в прак¬
тических задачах частично объясняется центральной пре¬
дельной теоремой, которая утверждает, что это распреде¬
ление вполне естественно появляется в результате суммар¬
ного действия большого числа независимых случайных вели¬
чин. Точнее, пусть х\, Х2, •••, хп есть п взаимно независи¬
мых случайных величин с произвольными и, возможно, раз¬
личными функциями распределения. Пусть pi и <у\ — среднее
значение и дисперсия случайной величины Х{, i = 1,... ,п. Рас¬
смотрим сумму случайных величин
п
X — ^ ^ QjjXj,
г=1
32

где di — произвольные фиксированные постоянные. Тогда сред¬
нее значение рх и дисперсия о\ случайной величины х имеют
вид
п	п
Их = ^2 афг,	.’ = £ afaf.
i=1	i=1
Последнее равенство справедливо в силу взаимной независимо¬
сти Х{	и	Xj	при	i	у£ j.	Центральная предельная	теорема
утверждает,	что при достаточно общих условиях распределе¬
ние суммарной случайной величины х при п —> оо стремится
к нормальному распределению с приведенными выше средним
рх и дисперсией
Свойства нормального распределения.
1. Если X и Y — два случайных вектора, совместная плот¬
ность вероятности которых имеет вид (3.39), то /(х), f(y)
и f(x/y) также имеют вид (3.39).
2. Если Y = L[X] — линейный оператор, то из нормальности
распределения вектора X следует нормальность распреде¬
ления вектора Y.
3. Для случайных величин, распределенных по нормально¬
му закону, понятия независимости и некоррелированности
совпадают.
4. Нормально распределенный вектор X с математическим
ожиданием рх и ковариацией Rx может быть получен с
помощью линейного преобразования
х — Sx • и рХ)
где и — нормально распределенный вектор с нулевым
средним значением и ковариацией, равной единичной мат¬
рице:
М [и] =0, М [и • иТ] = Е.
Для отыскания Sx надо решить задачу факторизации, т. е.
задачу представления Rx в виде
Rx = Sx • Sj при Rx > 0.
33

Матрицу S принято называть квадратным корнем из матрицы
Rx. Задача факторизации может быть решена неоднозначным
образом. Приведем два наиболее употребительных варианта.
1. Известно, что квадратичная форма ат Ra ортогональным
преобразованием /3 = Uа приводится к диагональному ви¬
ду:
aTRa = /3ТЛ/3, где Л = diag(Ab А2,..., Ап).
В нашем случае все А* > 0. Отсюда следует R = UAUT и
S = l/Лз, = diag (\/Ai, дЛг, •. •, V^) •
2. В качестве матрицы S выберем нижнюю или верхнетре¬
угольную матрицу вида
sh
0 ..
0
\
/
511
512
5ln
\
Sh =
«21
522
.. 0
SP
II
0
522
S2n
5nl
Sn2 • ■
. . Snn
/
\
0
0
Snn
/
Вычисление элементов матриц 5Н> 5в определяется
алгоритмами аналитического разбиения Холецкого [5].
Заметим, что треугольные матрицы широко используются
при решении задач оценивания.
Их очевидные свойства:
1. Сумма и произведение однотипных матриц является тре¬
угольной матрицей того же типа.
2. Обратная матрица от треугольной - есть треугольная мат¬
рица того же типа.
3.1.3.	Случайные процессы и их характеристики
Случайным процессом (случайной функцией) X(t,cj) будем на¬
зывать функцию двух аргументов: неслучайного параметра t и
элемента wGfi.
34

При фиксированном t величина X является случайным век¬
тором, при фиксации и; — детерминированной векторной функ¬
цией (реализацией случайного процесса). Параметр t далее бу¬
дем интерпретировать как время. Параметр t либо изменяется
непрерывным образом (процессы с непрерывным временем), ли¬
бо принимает дискретные значения, допускающие нумерацию
(процессы с дискретным временем).
Самый простой способ описать случайный процесс — это за¬
дать закон распределения процесса в каждый момент време¬
ни t. При таком способе вводится либо функция распределения
F(x,t), либо плотность вероятности f(x,t). В соответствии с
этим определяются математическое ожидание M[x(t)\ = fix(t)
О	О	~Г
и матрица ковариации M[X(t)X (£)] = Px(t,t).
Более полное описание случайного процесса получим, если
будем учитывать связь между значениями процесса в моменты
t\ = t и t2 = s, т.е. зададим совместный закон распределения
двух векторов X(t) и X(s):
F{xi,x2,t,s) = P{X(t) < х\,X(s) < х2} .
Для них можно определить матрицу
Px(t,s) = M[X(t),XT(s)],
которую будем далее называть матрицей корреляции, в отличие
от ее частного случая — матрицы ковариации Px(t,t).
Теория, основанная на описании процесса с помощью харак¬
теристик его вероятностного распределения для двух моментов
времени t и s, носит название корреляционной теории. В
дальнейшем все рассуждения будут проводиться в рамках этой
теории.
В приложениях важную роль играет принадлежность или
не принадлежность к трем классам случайных процессов. Эти
классы таковы:
• стационарные,
• марковские,
• нормальные (гауссовы) процессы.
35

В рамках корреляционной теории стационарными про¬
цессами называются такие, функции распределения которых
не зависят от выбора начала отсчета времени:
F(x(t),x(s)) = F(x(t + т), x(s + т)).
Для стационарного процесса
ixx — const, Dx = const, Px(t, s) = Рх{т), где т = t — s.
Компоненты вектора X в этом случае называются стационарно
связанными.
Процесс называется марковским, если вероятностные
свойства процесса в будущем полностью определяются вероят¬
ностными свойствами процесса в настоящем и не зависят от
поведения процесса в прошлом.
Пусть t\ < £2 < • • • < tn - п произвольных моментов време¬
ни. Тогда
P{X(tn) < Xn/X(tn-1) = xn_i,...,X(ti) = xi} ее
= P{X{tn) < Xn/X(tn-1)},
или, в терминах условной плотности,
/(%п/^п—1) • • • } 3^1) = /(з^п/з^п—l)*
Процесс называется нормальным или гауссовым, если
для любых моментов времени Ц, £2,- • -Дп соответствующие слу¬
чайные векторы имеют совместное нормальное распределение.
3.1.4.	Анализ случайных процессов
Цель последующего изложения — ввести понятия сходимости,
непрерывности, производной, интеграла применительно к слу¬
чайным функциям. Эти понятия, известные из области анализа,
нуждаются в существенном уточнении. Имеется три различных
определения сходимости последовательности случайных вели¬
чин.
36

Сходимость по вероятности.
Последовательность {АГп(о;)} сходится по вероятности к
Х(о;), если для любого е > О
Иш Р{\Хп-Х\ >е} = 0.
п—УОО
Сходимость в среднеквадратичном.
Последовательность {Хп} сходится в среднеквадратичном к
X, если
Иш М \\Хп -Х\2] =0.
п—УОО	1
Из сходимости в среднеквадратичном вытекает сходимость
по вероятности. Этот факт является непосредственным след¬
ствием неравенства Чебышёва:
Аналогично определяется непрерывность случайных функций.
Случайная функция x(t) непрерывна по вероятности в точ¬
ке £, если для Уе > 0
YimQP{\X(t + At) - X(t)| > e} = 0.
Случайная функция X(t) непрерывна в среднеквадратич¬
ном в точке t, если
дНто М [|X(i + At) - X(t)|2] = 0.
Оба определения непрерывности ничего не говорят о непре¬
рывности каждой реализации. Об этом говорит третье опре¬
деление сходимости, называемой сходимостью с вероятностью
единица или почти наверное. Но в дальнейшем нам это третье
определение не понадобится.
Теорема. Операции предела и математического ожидания
переставимы, если предел понимать в среднеквадратичном:
Иш М [Хп] = М\ lim Хп1 =М[Х].
п—У ОО	Ln—УОО	J
Теорема. Для того чтобы случайный процесс был непре¬
рывен в среднеквадратичном, необходимо и достаточно, чтобы
37

были непрерывны его математическое ожидание fj>x(t) и в диа¬
гональных точках функция корреляции Px(t,s) (при t = s).
Из непрерывности функции корреляции в диагональных
точках следует ее непрерывность всюду.
Производной в среднеквадратичном называется
предел в среднеквадратичном отношения
x(t + At) — x(t)
At
Определение означает, что
' x(t + At)
при At -» 0.
lim М
At—>-0
■x(t)
At
dxY
dt)
Теорема. Для того чтобы случайная функция была диф¬
ференцируема в среднеквадратичном, необходимо и достаточ¬
но, чтобы существовали производная математического ожида¬
ния и вторая смешанная производная функции корреляции в
диагональных точках.
Из существования второй смешанной производной в диаго¬
нальных точках следует существование вторых и первых про¬
изводных от этой функции всюду.
Поскольку производная является пределом в среднеквадра¬
тичном, операции дифференцирования и математического ожи¬
дания переставимы:
М
dx
= ^, м
'dS(t)
dx(s)
dt
dt ’
dt
dt
d2Px{t,s)
dtds
(3.40)
В частном случае стационарной функции получаем
dx	„ ,	d2Px(r)	п , ч
V = S'	Р»(м) =	—W~ = F'(T)'	т = ‘
s.
Интегралом f x(t)dt в среднеквадратичном называ-
а
ется предел в среднеквадратичном интегральной суммы
^	I
lim x(rk)At, At =	-,	(к	—	1)	•	At	<	т	<	к	•	At.
п—юо	п
38

Теорема. Для того чтобы процесс был интегрируем в сред¬
неквадратичном, необходимо и достаточно, чтобы существова¬
ли интегралы
ь ь
J fJ>x(t)dt,	J J Px(t, s)dtds.
а	а а
Очевидна перестановочность операций интегрирования и
математического ожидания.
3.1.5.	Процессы с ортогональными приращениями.
Белый шум
Процессы такого типа играют большую роль при спектральном
описании стационарных случайных процессов, при построении
стохастической модели динамических систем, в задачах анализа
и синтеза таких систем.
Введем процесс x{t) с нулевым средним M[x(t)\ = 0 опре¬
делив его следующими свойствами. Пусть £оДъ • • • Дь • • • ^тп ~
произвольные моменты времени и целое число га произвольно,
причем tk+1 > tk. Величина
Ax(tk,tk—i) = Ахк = x(tk) — x(tk—i)
составляет приращение процесса на интервале Дtk = tk — tk-1-
По определению, процесс будет процессом с независимы¬
ми приращениями, если Ах к не зависит от A Xj (к ф j) и для
У к Ах к не зависят от хо = x(to).
Если имеет место только некоррелированность, процесс на¬
зывается процессом с ортогональными приращениями.
Рассмотрим приращение процесса x(t) на интервале от t до
t + Ati
Ах = x(t + At) — x(t).
Выясним как связана ковариация этого приращения Рдх =
М [Да; • Да;1"] с ковариацией процесса Px(t) = М [x(t) • xT(t)\.
Имеем
РАх = М |^(a;(t + At) — x(t)) • (x(t + At) - z(£))Tj =
39

Px(t + At) + Px(t) - M [(x(t + At) - x(t)) • xT(t)] -
M |ж(£) • (x(t + At) — x(t))TJ
Ho
x(t + At) = x(t) + Ax,
и в силу определения процессов с ортогональными приращени¬
ями
М [(x(t) + At) - x(t)) • xT(t)] = Px(£),
M ^3c(t) • (x(t + At) - x(t))Tj = Px{t).
Отсюда
Pax = Px{t + At) - Px{t),
то есть ковариация приращения для процесса с ортогональны¬
ми приращениями равна приращению ковариации процесса.
В скалярном случае получим
DAx = Dx(t + At) - Dx(t).	(3.41)
Обозначим через R(t) производную по времени от ковариа¬
ционной матрицы процесса
R(t) = Px(t),
или, в скалярном случае,
dDx(t)
R(t)
dt
Последнее соотношение приближенно можно записать при до¬
статочно малом At
Daх ~ Rx(t)At,
т. е. дисперсия приращения процесса с ортогональными прира¬
щениями пропорциональна величине временного интервала At.
На первый взгляд это может показаться странным (если
не обращать внимания на условия дифференцируемости), по¬
скольку дисперсия — квадратичная характеристика. В самом
40

деле, рассмотрим процесс х(£), дифференцируемый в средне¬
квадратичном:
/ dx V
\Л)
At2.
Ах = x(t + At) — x(t) « xAt —> Ax & M
С другой стороны
DAx = M [Дя2] » At2M [х2] ,
то есть дисперсия приращения оказалась пропорциональной
квадрату приращения аргумента. Но процесс с ортогональными
приращениями, хотя и является непрерывным в среднеквадра¬
тичном, не имеет производной, поскольку
М [x(t + At) — x(t)] « Rx(t)At —> 0 при At -> О,
М
v(t -f- At) — x(t)
At
Rx(t)
At
—» oo при At —> 0.
Покажем, как можно интерпретировать величину dx, кото¬
рую нельзя понимать как обычный дифференциал в средне¬
квадратичном .
Предварительно напомним определение и свойства
6—функции. Эта функция широко используется при опи¬
сании линейных динамических систем и случайных процессов.
Она определяется следующими тремя условиями:
1.	5(х) = 0 при х ^ 0;
2.	J(x) = ос при х	= 0;
оо
3.	/	S(x)dx = 1.
—оо
Из условия 3 следует, что S—функция размерная величина.
А именно, ее размерность обратна размерности величины х.
В дальнейшем S—функции будут использованы для опи¬
сания функций корреляции Kx(t,s), обладающих свойством
Kx(ti$) = Kx(s,t), и	поэтому мы должны приписать
41

S—функции четвертое свойство - свойство четности: 5(х) =
6(-х).
Заметим, что четная 6—функция может быть определена
при помощи предельного перехода
1	_ х2
£(ж) = lim £*(ж, сг), где £*(ж, а) = - /— е .
л/27г а
Основное следствие из определения S—функции. Пусть <р(х)
гладкая функция, тогда
+оо	а+е
J (p(x)S(x — a)dx = J <р(х)5(х — a)dx =
—оо	а—е
где е как угодно малое положительное число.
Если 5(х) четная функция, то имеем
+оо	а+е
J <р(х)5(х — a)dx = J <р(х)6(х — a)dx =
— оо	а—е
а	а+е
J (f(x)S(x — a)dx + J ip(x)5(x — a)dx =
a—e
^(«) + \<p(a) = ¥>(«)•
Приращение процесса Аж на интервале (t, t + At) запишем
в виде
t+At
Аж
t
= J dx(u).
Тогда
t+At t+At
Pax = M [Ax • Джт] = J J M [dx(u)dxT(щ)] .	(3.42)
t t
С другой стороны,

Из сравнения (3.42) и (3.43) следует с учетом свойств
(5—функции
М [dx(u)dxT(иi)] = Rx(u)S(u — u\)dudu\.
Если ввести формальную производную,
получим
М [r(t) • rT(s)] = Rx{t)5(t - s).	(3.44)
Напомним, что в силу определения матрицы корреляции,
S—функцию в последних выражениях следует считать четной
функцией своего аргумента.
Сконструированный таким образом процесс r(t) называется
процессом типа белого шума или просто белым шумом. Ос¬
новное его свойство — некоррелированность во времени. Строго
говоря, такой процесс (3.44) — не более чем удобная абстрак¬
ция. Он не может быть реализован физически. Для его реали¬
зации требуется существование бесконечной скорости (абсолют¬
ной безынерционности, бесконечной энергии и т.п.).
3.2.	Стохастические модели линейных
динамических систем
Поведение динамических систем, которые описываются обык¬
новенными дифференциальными уравнениями, определяется
начальным значением всех фазовых координат (начальным зна¬
чением вектора состояния).
В стохастической ситуации задаются вероятностные ха¬
рактеристики начального состояния системы. Стохастическая
модель, описывающая состояние системы, должна позволять
предсказывать будущие вероятностные характеристики систе¬
мы. Для нас такими характеристиками являются математиче¬
ское ожидание и матрица корреляции вектора состояния. При
этом рассматриваются только линейные системы.
43

3.2.1.	Дискретный случай
Математической моделью системы служит вектор x(t), подчи¬
няющийся дискретному уравнению
Хк+1 = ФкХк + Чк,	(3.45)
М[чк\ — 0, M[qkq])=Qkhj, М[хкЧк\ = О, (Л-номер шага).
Здесь Skj - символ Кронекера.
Из (3.45) непосредственно следует уравнение для определе¬
ния математического ожидания = М [xk\:
№+1 = $кЦк-	(3.46)
Подставив (3.45) и (3.46) в выражение, определяющее ковариа¬
ционную матрицу Pk+1 = M[xk+iXk+iT], получаем ковариаци¬
онное (дисперсионное) уравнение:
Рк+1 = ФкРк$к + Qk•	(3.47)
о оТ
ХпХш
Легко показать, что матрица корреляции P(n, т) = М
определяется соотношением
Р(п,ш) = Ф(п,га)Рт.
Свойство. В стационарном случае (Q = const, Ф = const и
все корни характеристического уравнения |АЕ — Ф| = 0 по мо¬
дулю меньше единицы) при к —» оо можно говорить об устано¬
вившемся решении Роо дисперсионного уравнения (3.47). Оно
определяется из уравнения
Рх = ФРооФТ + Q.
3.2.2.	Непрерывный случай
Непрерывное стохастическое дифференциальное уравнение,
описывающее состояние динамической системы, будем рассмат¬
ривать как предел в среднеквадратичном следующего разност¬
ного уравнения:
Да; = А • xAt + Дг>,
44

Ах = x(t + At) — x(t), Av = v(t + At) — v(t),
где v(t) — процесс с ортогональными приращениями.
Обозначим производную ^ (Pv(t) = М [v(t)vT(t)]) че¬
рез Q. Имеем
М [Дг>Дг;т] = Q(t)At + o(At).
В пределе при At —> 0 получим
dx = A’xdt + dv,	M[dv(t)dvT (s)] = Q(t)dtds5(t — s).
(3.48)
Формально, если ввести производные ^ и можно записать
dx
— = Ах+ q,	(3.49)
где q — процесс типа белого шума:
М [g(t)] =0, М [g(i)gT(s)] = Q(t)S(t - s), M[x0qT(t)} - 0.
Уравнение для определения математического ожидания /хх =
М [х] непосредственно следует из (3.48) или из (3.49):
fix ~ A[J,Xi flx(tо) = flo-
Уравнение для определения ковариационной матрицы Px(t) =
М ^c(£)xT(t)J получим двумя способами.
1.	Примем за основу уравнение (3.48). Выпишем цепочку
простых соотношений:
Px(t) = М x(t)xT(t)j , Px(t+At) = M [x(t + At)xT(t + At)
x(t H- At) = x(t) + Ax(t),
Px(t + At) = M |x(t) + Дх(£) j |x(t) + Дж(£) j | =
= M |х(£)жт(£) И- M x(t)AxT(t)J + M j^\x(£)£T(£)j +
45

+м |^д£(г)д£т(£) .
Так как Ах = Ах At + Av, то
М |x(t)A:rTj = Px(t)AT(t)At, М |д£(£)ят| = A(t)Px(t)At,
М [д£(*)д£т] = A(t)Px(t)ATAt2 + QAt + o{At).
Таким образом,
Px(t “Ь At) =
= Px{t) + (APX + PxAT)At + QAt + APXAT At2 + o(At).
При At —У 0 из (1) получим
dPx
dt
'■ APX + PXA + Q.
(3.50)
Это уравнение называется дисперсионным или ковариа¬
ционным.
2. Примем за основу уравнение (3.49). Выпишем цепочку со¬
отношений:
^=Ax + q,	Px(t) = М [x(*)'cT(f)j ,
= М
dPx
dt
= АРХ + РХАТ + М
Поскольку
dx о -г о dx -г I
—хт +х— Ч =
dt	dt	1
ЭТ
q(t)x (t)
+ м
[£(*М*)Т] •
ТО
М
г
x(t) = $(t, <o)x(to) + У ф(*, r)q(r)dT,
to
t
x(t)g(t)TJ = J Ф(t,r)M [q(T)q(t)r]dT =
to
46

t
= J ^(t,r)Q{r)8(t — r)dr =
to
t
I
t—£ (e—>-0)
Отсюда следует уравнение (3.50).
Матрица корреляции Px(t,s) = М £(t)£T(s)j, как легко пока¬
зать, находится из соотношения
Px(t, s) = Ф(£, s)Px(s) (t > s).
В стационарном случае (Q = const, А = const и корни харак¬
теристического уравнения \ХЕ — А\ = 0 имеют отрицательные
действительные части) установившееся значение матрицы кова¬
риации может быть найдено из решения алгебраического урав¬
нения
АРХ + РХАт -j- Q = 0.
з.2.3. Дискретизация непрерывных случайных
процессов
Рассмотрим векторный случайный процесс x(t), описываемый
уравнением
х = A(t)x + q,	(3.51)
где х — n-мерный вектор состояния, q — n-мерный вектор —
процесс типа белого шума:
м [g(t)gT(s)] = Q(t)S(t - s).
Будем рассматривать процесс в дискретные моменты времени
t0, t\ = to At, ..., tfc = to + A;At, ...
и, соответственно, рассмотрим вопрос о получении дискретной
во времени модели системы, эквивалентной ее непрерывной мо¬
дели.
47

Определение. Дискретное и непрерывное во времени пред¬
ставления называются эквивалентными (в рамках корреляци¬
онной теории), если векторы состояния имеют одни и те же
математические ожидания и ковариационные матрицы в дис¬
кретные моменты времени.
Решая уравнение (3.51), получим
*fc+i
x(tk+1) = $(tk+i,tk)x(tk) + J $(tk+i,T)q(r)d,T,
где <l?(£fc+i,£fc) является переходной матрицей системы
Ф(£, t0) = A(t) Ф(£, tQ),	Ф(£0,	t0)	=	Е.
Введя обозначения
*к+1
Xk = x(tk), &k = ${tk+i,tk), Як= j ^{tk+uT)q{r)dT,
tk
получим дискретную во времени, эквивалентную модель
%к+1 = Фк%к “Ь Qki
где последовательность {^} — последовательность типа дис¬
кретного белого шума:
*fc+i
м [qkqj] = Q*kSki, Qk ~ J $(tk+i,T)Q{T)$T(tk+i,T)d,T.
tk
Замечание. При малом At обычно полагают
$(tk+i,tk) ^Е + ДtA(tk),
либо
At2
$(tk+1,tk) c?E + AtA(tk) + -y-A2(tk).
Тогда
Qt - Q(tk)£d.
48

Примеры.
1.	Процесс 1-го порядка.
Непрерывная модель.
х + \х = cry/2Xq(t),	М	[<q(t)qT(s)] = 1 • S(t — s),. (Л > 0).
Корреляционная функция Кх(т) процесса x(t) имеет вид
Zfc+1 = exp АЛ‘ xk + <г\Л ~ exp-2AAt qk, M [qkqi] = 1 • Ski■
2.	Процесс 2-го порядка.
Непрерывная модель.
х + 2ах + (Л2 + а2)х = 2<л/а(А2 + а2) q(t), (А > 0, а > 0),
[4(*)<?T(S)] = 1 *	—	s).
Корреляционная функция Кх(т) процесса x(t) имеет вид
Кх(т) = а2ехр Л,т|
Дискретная модель.
Кх{г) = а2 ехр а1т* (cos Ат + ^ sin А|т|^ .
Дискретная модель.
Ф =
(
м [qkqi] = Q*Ski, Q* = а2 [Е- Ф • Фт] .
49

Замечание. Возникающая задача моделирования дискретного
белого шума qk с заданной матрицей интенсивности Q* раз¬
решается, например, следующим образом. Факторизуем мат¬
рицу Q*
где Sq , для определенности, верхнетреугольная матрица. То¬
гда случайный вектор qk
будет иметь заданную ковариационную матрицу Q*.
3.3.	Алгоритмизация задачи оценивания
Последующее изложение посвящено задачам стохастически оп¬
тимального оценивания некоторой скалярной или векторной
случайной величины х с помощью скалярной или векторной ин¬
формации z. Оценка х реализует критерий, минимизирующий
в том или ином смысле ошибку оценки Ах = х — х. При этом
критерий отражает вероятностные связи, которые существуют
между измерениями z и оцениваемой величиной х.
Оценка может быть разовой (локальной) или меняющейся
со временем. Соответственно условно можно говорить о стати¬
ческих и динамических задачах оценивания.
Существует много различных методов стохастически опти¬
мального оценивания. Источники, излагающие соответствую¬
щую теорию с той или иной степенью детализации, весьма мно¬
гочисленны. Ниже мы остановимся на задачах, связанных с
тем, что называют обычно теорией калмановской фильтрации,
и примыкающих к ним.
Для суждения о качестве оценки обычно используются две
характеристики: несмещенность и состоятельность.
Оценка х называется несмещенной, если в среднем по ве¬
роятности она равна оцениваемой величине:
<3* = sqsг,
где
50
М[х] = х.

Пусть х — оценка, использующая результаты п измерений.
Оценка х называется состоятельной, если она сходится по ве¬
роятности к оцениваемой величине:
Р{\хп — х\ >£:}-> 0 при п —» оо.
В некотором смысле состоятельность является стохастическим
аналогом асимптотической устойчивости.
Отметим, что практически проще оказывается проверка бо¬
лее сильной сходимости — в среднеквадратичном.
Количественной характеристикой результата оценивания
служит оптимальность. Оценка оптимальна, если она наилуч¬
шая из всех возможных оценок с точки зрения некоторого за¬
данного критерия.
В связи с выбором критерия сделаем одно важное для при¬
кладных задач замечание. При построении алгоритма оценива¬
ния должны быть прежде всего четко сформулированы модели
всех инструментальных погрешностей измерителей, доставляю¬
щих первичную информацию, и модели возмущений, действую¬
щих на динамическую систему. Инженеры называют такие мо¬
дели обычно техническими условиями. Именно выбранная мо¬
дель погрешностей определяет критерий, а не наоборот.
Далее везде принято, что погрешности и возмущения име¬
ют нормальное распределение и характеристики такого распре¬
деления известны. Это приводит к квадратичному критерию
точности, который может использоваться в различных формах:
критерий минимума дисперсии ошибки оценки, критерий мак¬
симального правдоподобия, критерий ортогональности.
Ниже, в качестве основной формы выбран критерий ортого¬
нальности. Смысл его в рамках корреляционной теории состоит
в следующем. Оценка х, использующая некоторое измерение z,
должна быть такова, чтобы ошибка оценки Ах = х — х была не
коррелирована с исходной для оценивания информацией z:
М \AxzT 1 = 0.
Эвристически это означает, что вся содержащаяся в изме¬
рении z информация о величине х вошла в оценку х и, следо¬
вательно, в измерении z не содержится никакой информации о
неоцененной части величины х — ошибке Ах.
51

Ниже, в полном соответствии с гипотезой о нормальности
стохастических процессов, будут рассматриваться только ли¬
нейные алгоритмы оценивания.
3.4.	Решение переопределенных систем
линейных алгебраических уравнений.
Вероятностная интерпретация
Рассмотрим простейшую задачу решения системы линейных
алгебраических уравнений
z = Их,	(3.52)
где х — n-мерный вектор, подлежащий определению; z — из¬
вестный m-мерный вектор (вектор измерений); матрица Н име¬
ет размерность (га х п) и предполагается имеющей максималь¬
ный ранг (на практике имеет место чаще всего га>пи соответ¬
ственно система линейных уравнений (3.52) переопределена).
Рассмотрим два возможных подхода к решению этой зада¬
чи, первый связан с применением классического метода наи¬
меньших квадратов (МНК), второй основан на использовании
критерия максимального правдоподобия.
3.4.1.	Метод наименьших квадратов
Очевидный подход к решению задачи (3.52) таков. Поскольку
каждому вектору х Е Rn соответствует вектор невязки г Е R771:
г — z — Их,
то в качестве решения задачи (3.52) естественно выбрать такой
вектор х, который доставляет минимум (длины) нормы вектора
невязки г.
Тем самым приходим к классической постановке метода наи¬
меньших квадратов:
х = argmin J(x), J(x) = (z — Hx)T (z — Их) — rTr = ||r||2.
(3.53)
52

из которого следует
(3.54)
Матрица
Н+ = (ЯТЯ)-1ЯТ
называется псевдообратной. (Матрица ЯТЯ не вырождена, по-
скольку матрица Н максимального ранга.)
Вычисление псевдообратной матрицы Н+ в методе наимень¬
ших квадратов допускает рекуррентную процедуру, позволяю¬
щую уточнять решение (оценку) х от измерения к измерению.
Эта процедура будет рассмотрена ниже в разделе, посвященно¬
му дискретному фильтру Калмана.
3.4.2.	Критерий максимального правдоподобия
Придадим задаче решения системы линейных алгебраических
уравнений (3.52) вероятностный смысл. Будем считать, что
имеется вектор измерения г, такой, что
где г — случайный вектор погрешностей измерения, распре¬
деленный по нормальному закону с нулевым математическим
ожиданием и известной ковариацией R:
Замечание. Форма представления вектора измерения z (3.55J
имеет строгий смысл в отличие от записи (3.52,), в которой
знак равенства условен.
Функция плотности вероятности /(г) вектора г имеет вид
г = Hx + г,
(3.55)
м [г] = О, R = M [rrT] .
(3.56)
/(г) = С ■ exp [ ~^rTR V],
(3.57)
53

где
С = 2тг" ? (detR)~* .
Естественный подход к решению задачи (3.56) в рамках при¬
нятой модели помехи г таков. Поскольку имеет место (3.52),
то для любого значения вектора х случайный вектор z рас¬
пределен по нормальному закону с известной ковариацией R
и неизвестным математическим ожиданием Соответственно
функция распределения /(z,/xz) вектора 2 имеет вид
f(z, Hz) = с ■ exp -^(z-R~1(z-цг) .	(3.58)
Согласно критерию максимального правдоподобия, в каче¬
стве оценки Дг неизвестного параметра распределения (3.58)
принимается такое значение /хг, при котором функция плотно¬
сти распределения /(z, fiz) для полученного в результате опыта
(измерения) значения случайного вектора z достигает макси¬
мума.
Условие таx/(z,/xz), очевидно, эквивалентно условию
min J, где
J = (z — Нх)Т R-1 (z - Нх).	(3.59)
Тем самым, опять приходим к стандартной форме метода наи¬
меньших квадратов (см. (3.53)), где обратная ковариационная
матрица R-1 играет роль весовой.
Минимизация функционала J (3.59) приводит к следующе¬
му решению (см. п. (3.4.1)):
х = (HTR~1H)~1 HTRrlz.	(3.60)
При этом
Дх = х — х = х — (ЯтД-гЯ)-1 ЯТЯ-1 (Нх + г) = —Кг,
Pax = М [ДхДхт] = KRKr = (ЯТЯ"1Я)~1.	(3.61)
Очевидно, что при равноточных независимых измерениях
(R = сг2Е) из (3.60) вытекает знакомый результат (3.54).
54

Замечание. Пусть помимо измерений z привлекается точ¬
ная информация
z* = Gx, dimz*=l, I < п.	(3.62)
Существует 3 способа ее использования:
1. выразив из (3.62) I составляющих вектора х через
остальные, подставим их в уравнение (3.52) и тем са¬
мым уменьшим на I размерность вектора, определяемого
по МНК;
2. можно воспользоваться известным аппаратом множи¬
телей Лагранжа;
3. способ, целесообразный при практических расчетах: обра¬
зуем вектор w = aez*, где эе достаточно большое число.
Имеем
w = aeGx.	(3.63)
Добавим уравнение (3.63) к уравнению z = Нх. Оцен¬
ку вектора х будем искать как решение совокупной си¬
стемы, в которой вектор измерения имеет размерность
т + 1. Это прием эквивалентен тому, что компонентам
вектора z* приписываются малые погрешности со сред¬
неквадратичными отклонениями
3.4.3.	Задача сглаживания экспериментальных данных
методом МНК при помощи кубических сплайнов
В качестве примера задачи, приводящей к решению переопреде¬
ленной системы линейных алгебраических уравнений, рассмот¬
рим задачу сглаживания экспериментальных данных при помо¬
щи кубических сплайнов.
Пусть g(t) функция времени, измеряемая в дискретные мо¬
менты tfc:
^0>	•	•	•	>	tk-\-1)	• • • »	=	tk-\-1	tfc = COTlst.
Результат измерения
Zk =g{tk) + rk,
55

где rk - погрешность измерения.
Общее число измерений Zk функции д в моменты времени
{tk} равно N + 1. Требуется построить оценку g(t) на всем от¬
резке [to,*J\r]-
Отрезок [to, tjv] разбивается п узлами Т*:
При этом число п узлов разбиения отрезка [to, t^v] меньше
чем число измерений: N + 1 >77. Обычно для практических
задач N + 1 >> п.
Узлы выбираются так, чтобы им отвечали какие-то изме¬
рения. Измерения в узловой точке Ti удобнее относить к преды¬
дущему отрезку [Ti-i—Ti\. Число измерений на каждом участке
кроме первого обозначим через т. Тогда первый участок [Ti, Т2]
будет содержать т+ 1 измерений. Имеет место соотношение
771 + 1 + (п — 2)777 = N + 1 или 77(777 — 1) = N.
Введем безразмерный интервал между измерениями:
На каждом участке [Т^_ 1 —7*] введем безразмерное время т:
Основной элемент сглаживания при помощи кубических
сплайнов - определение базовых функций. Вводятся 4-е базо¬
вые функции:
так, что
= to, • • •, Тп = ijy,	=	7i_|_i	-Тг	=	const.

Линейная независимость функций v?i(r), <р2(т),	<Р4(т)	сле~
дует, например, из однозначности определения коэффициентов
С1,С2,сз,С4 в соотношении
(т) + Ci+2(p3(r) •+- Сг+з</?4(т), t G [Ti,Ti+i],
Непосредственной проверкой легко убедиться в том, что g(t)
непрерывна вместе со своей первой и второй производной в уз¬
ловых точках и, стало быть, всюду.
Для определения коэффициентов
имеем следующую переопределенную систему, имеющую так
называемую ленточную структуру.
Интервал [Ti,T2]:
zo	=	cnpi(0) + С2У>2(0) + сз<£з(0) + С4У>4(0),
Zi	=	ci(fi(Ar) 4- c2ip2(Ar) 4- с3(/?з(Дт) + с4<^4(Ат),
Z2	=	cnpi(2Ar) + с2(р2(2Ат) 4- с3^3(2Дт) 4- с4^4(2Дт),
Zm	=	Ci(pi(l) 4- C2(f2(l) + Сз^з(1) + C4^4(l).
Zm+l = c2<£i{At) 4- С3<Р2{Ат) + Сз^з(Аг) 4- същ{Ат),
Z2m = C2<Pl(l) + С3<р2( 1) 4- С4у?з(1) + С5^4(1).
Интервал [Tn_i,Tn]:
Ciy?i(r) 4- с2(р2(т) 4- с3^з(т) 4- c±ip±(r) =
ai 4- агт 4- а3т2 + а4Т3.
Оценка д функции д ищется в виде
57

Zn = C„_iV?i(l) + Сп¥>2(1) +с„+1</?з(1) -\-сп+2<Ра(1).
Решение этой системы обычно находят по методу наимень¬
ших квадратов.
Заметим, что алгоритм сглаживания может успешно при¬
меняться, если некоторые из измерений или даже целые блоки
таких измерений отсутствуют.
3.4.4.	Критерий ортогональности и критерий
условного среднего
Далее решается задача оценивания, которая служит основой
многочисленных модификаций метода наименьших квадратов.
Из нее, в частности, следует и предыдущий результат.
Рассматриваются два случайных вектора х(п х 1) и
z(m х 1) с известными характеристиками:
цх = М[х],	= M[z\,
Рхх = М [х£т] , Рг
М
о о”Г
ZZ
PXZ = M
о о~Г
XZ
Указанные характеристики будем называть априорной инфор¬
мацией.
Ставится задача определения (построения оценки) скаляр¬
ной величины а = стх, где с(п х 1) — известный вектор.
В рамках априорной информации в качестве оценки а ве¬
личины а естественно принять ее математическое ожидание,
т.е. положить а = ст//х. Ошибка оценки Да = а — а. Ме¬
рой ошибки оценки служит дисперсия D[Aa) = М
М [стххтс
‘1-
>Ti)1 -
С -PrrrC.
Пусть теперь измеряется вектор z и резуль¬
тат измерения	известен. Поскольку	известна	ме¬
ра	линейной	связи двух векторов	х и z	—
матрица Pxz, возникает задача получения оценки величи¬
ны а с учетом измерения z.
Оценку будем искать как линейную, несмещенную и опти¬
мальную при условии, что критерием оптимальности служит
минимум дисперсии ошибки оценки Да.
58

Линейность. Оценка а ищется в виде
а = 71T/ix + 72Т/хг + 73Tz,
где векторы 71/у2 и 73 подлежат определению.
Несмещенность. Условие несмещенности означает, что в
среднем оценивание приводит к точному результату, т.е.
М[Аа] = 0. Имеем
М[Да] = M[cTx-71T/Xx+72TMz+73T^] =
=	(с-	7i)%x-(7*+73)V
Отсюда 71 = с, 72 = —73. Переобозначим 73 = х, тогда
а - стМх + x(z — fiz).
Вектор ус подлежит дальнейшему определению.
Оптимальность.
£>[Да] = М j^(Aa)2j = сТРххс + xTPzzx - 2cTPxzx.
Условие оптимальности имеет вид
dD[Aa]
дк
откуда следует
= 0,	(3.64)
Хт = стРХ2Р-1 = ст К,	(3.65)
где К = PxzPzzl-
Если в качестве оценки х вектора х выбрать
X = fix + K{Z - fJLz),	(3.66)
то при всяком с оценка а = стх оптимальна.
Мерой ошибки оценки Ах = х — х вектора х служит кова¬
риация Рхх = М [ДхДят] = М (^х — Kz^j — Kz^j j . После
выкладок с учетом (3.65) получаем

Мерой ошибки оценки величины а служит дисперсия
D[Aa] = стРххс.
Величины х и Рхх можно интерпретировать как апостери¬
орные характеристики (условные, уточненные с учетом изме¬
рения z). Результаты, сведенные в формулы (3.65), (3.66), поз¬
воляют сформулировать критерий оптимальности (3.64) в виде
критерия ортогональности: оценка х оптимальна, если ошибка
оценки ортогональна (некоррелирована) измерениям:
В самом деле, из (3.68) следует (3.65) с учетом того, что Ах =
Как уже говорилось, критерий (3.68) служит отражени¬
ем простого эвристического рассуждения: оптимальная оценка
должна вобрать в себя всю информацию об оцениваемой вели¬
чине, содержащуюся в измерениях, и, следовательно, ошибка
оценки не должна быть связана с измерениями.
Ввиду важности полученного результата дадим ему несколь¬
ко иную интерпретацию, напрямую используя нормальность за¬
кона распределения случайных векторов х и z. Оценка х будет
определяться как условное математическое ожидание.
Пусть случайный вектор
распределен по нормальному закону с характеристиками
Функция плотности вероятности f(u) = f(x,z) вектора и
записывается в следующем виде:
М Дх5т] =0.
(3.68)
(3.69)
60

Пусть вектор z доступен измерению. Требуется уточнить
априорную характеристику М [х] случайного вектора х и опре¬
делить соответствующую ковариацию.
Для решения этой задачи естественно воспользоваться фор¬
мулой условной вероятности
Л./.) -
где функция распределения f(x,z) определяется соотношения¬
ми (3.70), а функция f(z) имеет вид
/и=“р В - фТр-ф - В •
Обозначим через fix/z и Px/Z условное математическое ожи¬
дание и условную ковариационную матрицу. Тогда
Пф) =	“Р{'5 (* '^‘)Тр^{х- ''*/*)}'
Рх/г = М[(Х- Цх/2) (X - flx/z)T jz .
Далее нам понадобится соотношение
р-1 _ ( Рхх Рхг \ 1 _ N ( Nxx Nxz N
V Pzx Pzz J	\	Nzx	Nzz J	’
где
NXX = (Pxx - PxzPzzPzx) 1 , Nzz = {Pzz - PzxPxxPxz) 1 ,
Nxz = - {Pxx - PxzPzzPzx) 1 PxzPzz-
Справедливость последних соотношений непосредственно сле¬
дует из равенства PUN = Е.
Применяя формулу Байеса и приравнивая коэффициенты
при одинаковых степенях х, получаем уже известный результат
ftx/z = Мх + PxzPzz (Z ~~ №z) j
61

Px/z — Pxx PVzPZZ Pzx — Pxx PXzPzz PXZ'	(3-71)
Рассмотрим частный случай, когда вектор z линейно зави¬
сит от вектора я, т.е.
z = Нх + г,
Н — известная матрица, г — случайный вектор с нулевым сред¬
ним, трактуемый как погрешность измерения.
Будем считать заданными априорные характеристики
Их = М[х], Рхх = М |ххт] , R = М [ггт] .
Кроме того, предположим, что вектор г некоррелирован с оце¬
ниваемым вектором х:
М [хгт] = 0.
Целесообразно ввести обозначения, в которых априорность
фиксируется верхним индексом «—>, апостериорность — верх¬
ним индексом «+»:
2/	—	Prg	—	Рхх'
После очевидных выкладок имеем
Р- =PZZ = HP~HT + R, Pxz = Р~НТ, !MZ = Нцх = Нх~.
Взамен соотношений (3.65), (3.66), (3.67) получим
х+ = хГ + К (z — Нх~),
К = Р~НТ (НР~НТ + R)-1,
Р+ = Р--КНР-.	(3.72)
3.5.	Дискретный фильтр Калмана
Дальнейшее изложение посвящено динамическим задачам оце¬
нивания, когда поведение во времени оцениваемой величины х
описывается линейными стохастическими (непрерывными или
дискретными) уравнениями. Будет использоваться только одна
форма критерия оптимальности — критерий ортогональности.
62

Обсуждение указанного круга задач начнем с дискретного слу¬
чая.
Термин «фильтрация» означает, что алгоритм может быть
интерпретирован как алгоритм выделения полезной информа¬
ции на фоне помех.
3.5.1.	Алгоритм дискретного фильтра Калмана
Задача в ее дискретном варианте ставится следующим образом.
Поведение динамической системы подчиняется уравнению
Xj-|-i = Ф j Xj -j- Qj,
где Xj — вектор состояния в момент tj; qj — дискретный белый
шум с известной интенсивностью — ковариационной матрицей
Qj; Фj — известная переходная матрица.
Заданы = М [т0], Ро = М ^(хо —	) (жо ~	)Т]	• ® м0_
менты времени j = 0,1,2,... поступает информация zo, zi,...:
Zj = HjXj + rj.
Погрешность информации rj — некоррелированный во време¬
ни вектор с нулевым средним и известной интенсивностью Rj\
М	= RjSij. Кроме того, полагаем М	=	0.
Требуется в каждый момент j получить оценку Xj, удовле¬
творяющую условию линейности, несмещенности, оптимально¬
сти. Введем обозначения:
xj — оценка вектора х в момент времени j, исполь¬
зующая измерения zo,...,Zj-i (априорная оценка х~ =
M[x/z0,z
х* — оценка вектора х в момент времени j, использу¬
ющая измерения zq, ..., Zj-i, zj (апостериорная оценка xj =
M[x/z
Р~ — ковариация ошибки априорной оценки вектора х в
момент j;
р7 — ковариация ошибки апостериорной оценки вектора х
в момент j.
Из предыдущих пунктов непосредственно следует, что про¬
цедуру оценивания можно представить в виде повторяющихся
63

от момента к моменту циклов, каждый из которых состоит из
двух этапов — этапа коррекции и этапа прогноза.
Этап коррекции — переход от априорных оценок к апо¬
стериорным в момент j. В соответствии с формулами (3.72)
= xj + Kj (Zj - Hjxj),
Kj = P~Hj {HjPfHj + Rj)1,
P+ = {Е-К^)Р~.	(3.73)
Этап прогноза — переход от апостериорных оценок х+ в
момент j к априорным оценкам х~ в момент j + 1. Имеем
*7+1 = ф**£,
Р~+1=Ф^+фТ +Qj.	(3.74)
Кроме того, для включения начального цикла надо положить
х0 = ^о> Р0 = ро-	(3.75)
Алгоритм оценивания, описываемый соотношениями (3.73)-
(3.75), носит название дискретного фильтра Калмана.
Пояснение. Из самой процедуры получения оценок сле¬
дует, что ошибка оценки Ах^ = Xj —х^ ортогональна текущему
измерению Zj\
М
Но она также ортогональна всем предыдущим измерениям.
Проиллюстрируем этот факт на примере двух изме¬
рений. Оценка х$ в начальный момент строится из условия
= 0, где Дх^ = х0-xj, 20 = z0 - ЯоМх,
рьх - априорное значение математического ожидания величины
XQ.
Если
xi = Ф0х0 + q0, zi = Hixi + n,
64
M
л + °Т
Ax£z0
Дх+£
— 0, где Zj = Zj — HjX-

то оценка х^ в момент времени t\ строится по условию
М
Ах+гг
= 0,
где
Axf
■ Х\ Х‘
Z-+
Z\ — Z\ — Н\хх .
Покажем, что величина Ахх ортогональна измерению zq. Для
этого ошибка Дя]1" должна иметь следующую структуру:
Д^1~ = AiAx^ 4- A2qo + A3ri.
Из выше написанного следует:
Axf = х\ —х\ — х\ —х\ — K\Z\— H\Xi ,
Xi-Xi = Ф0х0 + q0 - Ф0х£ = Ф0Ах£ + q0,
Z\ — HiXi = H\X\ -f ri — Hix^ = H\Ax~[ 4- ri =
= #1 (ФоДх^ + q0) + ri.
Таким образом мы показали, что величина Дх* имеет требуе-
мую структуру. Отсюда по индукции следует
Дж+ _L Zj-i _L Zj-2 ± ... ±zq.
Возможно обобщение дискретного фильтра Калмана на слу¬
чай коррелированности шумов системы и шумов измерений.
Пусть дополнительно М \qjrj] = Lj, тогда соотношения, опи¬
сывающие этап коррекции остаются прежними, а ковариацион¬
ное соотношение этапа прогноза будет таким
Р~+1 = Ф,Р+Ф] + Qj - bjKjL] - LjKjФ].	(3.76)
В частном случае, нередко встречающемся на практике, при
Qj = 4j + ®jrj, М [q*rj] = 0,
получим

- Ф,-ед-фТ - Ф,.Дл-к7ф7,	(3.77)
где Q* = М [q*q*T].
Некоторые свойства дискретного фильтра Калмана.
Стационарный случай:	= Ф = const, Qj = Q = const,
Rj — R — const.
Поскольку Q > 0, то существует представление (неоднознач¬
ное) Q = BBT, где матрица В называется квадратным корнем
из матрицы Q.
Теорема. Пусть пара (Ф, Н) наблюдаема, пара (Ф, В) управ¬
ляема. Тогда при бесконечном времени наблюдения:
1) существуют Р^, Р+, Коо, определяемые соотношениями
Р- = ФР+ФТ + Q,
Р+ = Р--КооНР-,
Коо = P-HT{HP~HT+R)
2) уравнения ошибок относительно величин Дxj = Xj — х~,
Ах* = Xj - Xj’
Ax~+1	= Ф Axj+qj,
Axf = (E — KooH) Ax~ — KooTj,
таковы, что при qj =0, rj =0, выполнено
AxJ, Axj~ —> 0,	{j -)• oo).
Более подробно случай стационарности будет рассмотрен в
разделе, посвященном непрерывному фильтру Калмана.
Пример. Пусть xj+i = Xj + Zj = Xj + 7^, Q = 1, R= |.
p-
Здесь Ф = 1, H = 1, Kj =	-	Найдем	установившееся
значение К и вид фильтра:
Р- = 3-	=	1	к	=2-
ОО 2’	ОО	2’	00
.+	.+ __	2 , __ч
хз+1 = ХТ ’	ХТ = хо +	з W ~хо)-
66

Уравнение ошибок таково:
1 л	2
Axj+i = -Axj + qj - -rjy
при этом однородное уравнение
Axj+l = \Axj
асимптотически устойчиво.
Полученные результаты, связанные с дискретным фильтром
Калмана, позволяют придать иную форму рассмотренному ра¬
нее методу наименьших квадратов в задаче решения переопре¬
деленной системы алгебраических уравнений z = Нх.
Рассматривая компоненты Zj как следующее одно за дру¬
гим измерения, методу наименьших квадратов можно придать
рекуррентную форму, интерпретируя индекс j как дискретное
время.
Будем считать, что в начальный момент .7 = 0 априорной
оценкой служит величина хо = 0. В качестве начальной кова¬
риации вектора х примем величину Ро = где ае достаточно
большое число (ае 1).
Далее используем процедуру фильтра Калмана, исключив
из нее этап прогноза. Будем считать, что все измерения некор-
релированы между собой и стохастически равноточны. Поло¬
жим для определенности, что М [r|] = 1.
Первый шаг:
х1 = х° + К\ (zi — hjx°) ,
Рг = Ро - PoM^a/ii +1)_1
Ki = P0hi (hJPohi + iy1.
Второй шаг:
x2 = x1 + K2 (z2 — h^x1),
P2 = P1-P1h2(f^P1h2 + l)~1h^P1,
K2 = P\h2 {h2P\h2 + l)-1.
67

j-ый шаг:
rjj	—	%3~1 _j_ J£.	_ foi x-7-1) ,
Pj	=	Pj—l — Pj—lhj {hj Pj—lhj + l)	hj Pj~i,
ifj	=	Pj-ihj (hJPj-ihj + l) 1.
Последний шаг j = т. Очевидно, что при достаточной боль¬
шой величине ае, хт с высокой степенью точности совпадает с
оценкой х = (ЯТЯ) 1 НТ z.
3.6.	Реализация дискретного фильтра Калмана
методом квадратного корня
На практике обычно используется численная модификация дис¬
кретного фильтра Калмана, основанная на так называемом ме¬
тоде квадратного корня.
Метод квадратного корня состоит в том, что ковариацион¬
ная матрица Р представляется либо в виде разложения Р =
55т, где S - верхне- или нижнетреугольная матрица, либо в
виде Р = UDUT, где U - верхне- или нижнетреугольная мат¬
рица с единицами на главной диагонали, D = diag (di, d2,
причем dj > 0. При таком представлении, каковы бы ни были
невырожденные матрицы S и D, матрица Р будет положитель¬
но определенной.
Метод корня в частности позволяет избежать потери по¬
ложительной определенности ковариационной матрицы Р на
этапе коррекции при вычитании, когда разности соизмеримы
с квантом дискретизации по уровню. В соответствии со сказан¬
ным, всюду в соотношениях, описывающих фильтр Калмана,
заменим ковариацию Р на произведение SST, где S = Р1/2
треугольная матрица.
Рассмотрим сначала случай скалярного измерения:
Zj = hjxj + т-j, M [rj] = 0, M [rirj] = RjSij.
Этап коррекции (нижние индексы пока опускаем). Имеем
Р+ = Р~ - P~h (hTP~h + R)~1 hTP~,
68

Х+ = Х~ + К (z - hTx~), К = p-h(hTP~h + Л)-1.
Положим Р~ = SI_5_T, Р+ = S+S+T. Получим
S+S+T = S~S~T - S~S~Th (hTS~S~Th + R)1 hTS~S~T,
К = S~S~Th (hTS~S-Th + Л)-1.
Введем обозначения f = S~Th, a = fTf + a2, a2 = R. Оконча-
тельно получим
1
5+ = S~
E--ff
a
T
1/2 1
К = —S~f.	(3.78)
a
Если z - вектор и матрица R - диагональна (составляющие
вектора г некоррелированы между собой), этап коррекции за¬
ключается в последовательной поскалярной обработке каждо¬
го измерения, причем апостериорное значение S каждого шага
служит априорной информацией следующего шага. Результат
обработки последней координаты вектора z служит исходной
информацией для перехода к этапу прогноза.
Замечание. Если матрица R не диагональна, то невырож¬
денным преобразованием С вводится вектор
z* = Cz = СНх + Сг = Н*х + г*
такой, что координаты вектора г* будут некоррелированы.
Например, матрицу С можно выбрать так, чтобы
M[r*r*T] = Е. Имеем в этом случае M[r*r*T] = CRCT = Е.
Откуда
С = R~1/2
— обратная матрица квадратного корня R1^2.
Этап прогноза.
Xj+i =	>	Pj+i	—	4-	Qj,
Sb.1s£1-*JSfS}T*J +Q,.
69
или

Обозначим Wj = $jSj~. Тогда
S-+1 = [WjWj + Qj]1/2 •
Модификация метода корня при помощи U-D разло¬
жения.
Положим Р = UDUT, где D - диагональная матрица с по¬
ложительными диагональными элементами, U - верхне- или
нижнетреугольная матрица с единицами на главной диагонали.
Структуры возникающих при таком представлении алгоритмов
во многом идентичны описанным выше структурам.
Этап коррекции (случай скалярного измерения).
Р~ = U~ D~U~T, Р+ = U+D+U+T.
Получим
х+ = х~ + U~vAz/a, U+ = U~U, D+ = D,
(г>- -	=	UDUT,	Az	=	z-hTx~,
f = U~Th, v = D~f, a = fTD~f + <rl
В случае коррелированных векторных измерений z использует¬
ся прием предыдущей модификации метода корня.
Этап прогноза.
Х~ = фх+, U~D~U~T = Фи+Р+и+тФт + Q.
3.6.1.	Численная реализация метода квадратного корня
фильтра Калмана для S-верхнетреугольной
формы
Опишем алгоритмы численной реализации метода корня для
его 5-верхнетреугольной формы [5].
1.	Алгоритм извлечения квадратного корня из матри¬
цы.
Р = SST, S = {Р}1/2 ,
70

S - верхнетреугольная матрица.
Обозначим элементы матриц Р(пхп) и S(nxn), соответственно,
через Pij и Sij.
1.	Вычисляем
2.	Численный алгоритм метода корня для этапа обра¬
ботки измерений. Поскалярная обработка.
Основу численной реализации этапа обработки измерений
(см. (3.78)), составляет возможность явного представления кор¬
ня
фигурирующего в (3.78), в виде функции от векторного пара¬
метра / и скаляра а. В конечном итоге это приводит к ком¬
пактному рекуррентному алгоритму.
В случае верхнетреугольной S'-формы метода корня он име¬
ет следующий вид.
2.	Для г = 1,2,... п — 1 вычисляем
Pin
&пп
3.	Для j = п — 1, п — 2,..., 2 вычисляем
О
г > j
п
Pij— 52 SikSjk
k=j +1
4.	Вычисляем
n
71

Введем обозначения
/ = S~Th =
(h \
/2
е0 -
О \
V 0 /
do - а2,
\ /» /
5- = (SfS2- •••5-),	5+ = (5+5+ ...5+). (3.79)
Здесь 5^ , 5^" - fc-ые столбцы матриц 5“ и S+, ео - нулевой
вектор размерности п.
Для к = 1,2,..., п вычисляем
dk
dfc-i + /fc> Ьк—
dk-1
dk	\J	dk-idk
&k — &k—i “Ь $k fk't ^~k ~ bk e^—iC/j.
После n итераций искомый квадратный корень имеет вид
s+ = (s+s+...s+).
Апостериорная оценка х+ вектора состояния х определяется
fk
выражением
х+ = х + е
Az
п~Т~
dn
= х + К {z — hTх ) ,
где
К
dn
калмановский коэффициент усиления.
3.	Численный алгоритм метода корня для этапа прогноза.
Для общности рассмотрим модификацию задачи дискретной
калмановской фильтрации, когда рассматривается следующая
модель динамической системы:
X j _|_ \ — Ф j X j "I- Bj Qj,
где Bj - известная матрица - матричный множитель при невы¬
рожденном шуме qj:
м = Qjfijki Qj > о»
72

Тогда дисперсионное уравнение этапа прогноза фильтра Кал-
мана примет вид:
P-+1 = Ф3Р+Ф] + BjQjBj.	(3.80)
Отправной точкой при разработке численной реализации со¬
отношений (3.80) для S-формы метода корня, является факто¬
ризованное представление этого уравнения (индекс j опускаем):
р- = S~S~T = Ф5+5+Тфт + В SqSjBT = ААТ,
где
А = ( Ф5+ BSq ),
a Sq - квадратный корень из матрицы Q:
Q = Sqsj.
Главная идея соответствующих вычислительных алгорит¬
мов заключается в определении такой ортогональной матрицы
Т((п + пя) х (n + nq)) (nq - размерность вектора q), которая
приводит коагулированную матрицу А к блочно-треугольному
виду
АТ= (5*: О).
Здесь S*(n х п) - треугольная матрица такого же типа, что и
матрица S+, О - нулевая матрица размерности (n х nq). Оче¬
видно, что матрица S* и является искомым квадратным корнем
s-.
Матрица Т, в частности, может быть определена с помо¬
щью процедуры ортогонализации Грама-Шмидта, примененной
к строкам матрицы А.
Компактная запись соответствующего вычислительного ал¬
горитма имеет следующий вид. Введем обозначения:
Ат = (*S+’:BSq)T = (jiiAa-.-A,),
*5”	= КЛ’ SqS* = Q.	(3.81)
73

Здесь Ak - к-ый столбец матрицы Ат размерности n + nq.
Для к = п, п — 1,..., 1 вычисляем
= \jAlAk,
Ак
Щ = —,
s
fcfc
А?
Здесь «— означает замену старого значения j-го столбца Aj
матрицы А новым значением.
3.7.	Непрерывный фильтр Калмана
Рассматривается непрерывная линейная динамическая систе¬
ма, поведение которой описывается уравнением
x(t) = A(t)x(t) + q(t).	(3.82)
Возмущения q(t) представляют собой нормальный случайный
процесс типа белого шума:
М [g(t)] = 0, М [g(i)gT(s)] = Q(t)6(t - s), Q(t) > 0.
Заданы также
x(to) = М [ж(£0)]>	P(t0)	=	М [(x(*o) - x(t0)) (x(t0) - х(£0))Т] > 0.
Возмущение q(t) предполагается независимым от начального
состояния системы x(to):
М [x(to)qT (t)] =0	Ш > to.
Пусть информация
z(t) = H(t)x(t) + r(t)	(3.83)
поступает непрерывно на интервале [to,t].
Погрешность информации r(t) — случайный нормальный
процесс типа белого шума:
М [г(^)гТ(5)] = R(t)5(t — s), R(t) > 0.
74

Кроме того,
М [x(to)rT(t)] =0 \ft > to, M [r(£)gT(s)] =0 Vt,s.
Требуется в любой момент времени t > to определить оценку
x(t), удовлетворяющую условиям линейности, несмещенности и
ортогональности.
Линейность. Оценку будем искать в виде
х = Tix + Г2z,
где матрицы Гi и Г2 подлежат выбору.
Несмещенность. Запишем уравнение ошибок оценки отно¬
сительно величины Ах = х — х:
Ах = Ti Да: + (А — Ti — Г2 Н)х + q — Г2Г.
Требование М [Да:] = 0 приводит к условию
A — Ti— Г2Я = 0 или	Т1=А-КН,	К = Т2.
Таким образом, уравнение для оценки х имеет вид
х = Ах + К (z — Нх).
Матрица К подлежит определению.
Ортогональность. Условие ортогональности имеет вид
М |^a:(t)!T(s)J =0, to < s < £,
где Ах служит решением уравнения ошибок
Ах = (А — КН) Ах + q — Кг	(3.84)
и
z(s) = z(s) — M[z(s)\ = z(s) — Нх = H(s)Ax + r(s).
Отсюда
Px(t, s)HT(s) + M [Ax(*)rT(s)] = 0.	(3.85)
75

Запишем решение уравнения (3.84):
t
Ax{t) = Ф(Мо)Ля(*о) + У $(t,T)[q(r) - К(т)г(т)]<1т,
*0
где Ф(£До) — переходная матрица соответствующей линейной
системы.
Далее получим
t
to
t
= —J Ф(£,т)КГ(т)Д(т)£(т — s)dr = — Ф(£, s)K(s)R(s).
to
Уравнение (3.85) теперь запишем так:
Px{t, s)HT(s) - Ф(£, s)K{s)R{s) = 0.
Устремляя s —> t, находим K(t):
K(t) = Px{t)HT(t)R~1{t).
Ковариационная матрица Px(t) = M [Ax(t)AxT(t)] находит¬
ся как решение дисперсионного уравнения, соответствующего
уравнению (3.84),
РХ = (А- КН) Рх + РХ{А - КН)Т + Q + KRKT. (3.86)
Подставив выражение для матричного коэффициента усиления
К в (3.86), получим
Рх = АРХ + РХАТ + Q- PxHTR~lHPx.	(3.87)
Уравнение вида (3.87) в теории дифференциальных уравнений
носит название уравнения Риккатпи.
Алгоритм оценивания, описываемый соотношениями
х = Ах + К (z — Нх),	x(t0) = хо,
76

K(t) = Px(t)FT(t)fl-1(<),
Px = APX + PXAT+ Q-PXHTR~1HPX,	Px{to) — Po>
(3.88)
носит название непрерывного фильтра Калмана.
Нетрудно привести обобщение фильтра Калмана для слу¬
чая, когда шумы в системе q(t) и шумы наблюдения r(t) корре-
лированы. Пусть имеет место условие
М [r(t)gT(s)] = W(t)S(t — s).
Повторяя предыдущие выкладки, получим
K(t) = [Px(t)HT + WT(t)] R-\t),
Px = APX + PXAT + Q-К [HPX + W].
Алгоритм, доставляющий оценку вектора состояния в соответ¬
ствии с формулами (3.88), вырождается при R = 0 (измерение
не содержит шума г). Такая ситуация возникает, например, ко¬
гда коррелированный шум в измерении представляется в виде
формирующего уравнения и включается в вектор состояния.
Один из удобных минженерных” вариантов решения задачи со¬
стоит в том, что вводится некоторая матрица Д, соответствую¬
щая фиктивному шуму г малой интенсивности:
М [ггТ] = еЕ {е > 0 — малый параметр).
Далее выясняются некоторые свойства алгоритма (3.88) и,
в частности, уравнения Риккати.
3.7.1.	Представление уравнения Риккати в виде
линейных уравнений большей размерности
Для краткости обозначим НтR~lH = D > 0. Прежде всего
заметим, что уравнение
Р = АР + РАТ + Q - PDP	(3.89)
при Q = 0 сводится к линейному путем перехода к обратной
матрице N = Р-1:
N = -NA - АтN + D.
77

Рассмотрим совокупность соотношений, определяющих матри¬
цы ФЬ \&2‘
Фх	=	A*i + Q*2,	®i(to) = JFb,
Ф2	=	DV1 - ЛТФ2,	*2(^0) = Е,	(3.90)
и определим матрицу Р* соотношением
фх = Р*Ф2.	(3.91)
Покажем, что соотношения (3.90)	эквивалентны уравнению
(3.89) в том смысле, что задают одну и ту же функцию P(t) =
P*(t). Продифференцируем соотношение (3.91):
ф1 = Р*ф2 + р*ф2
и подставив выражения для производных Фх и Ф2 из соотноше¬
ний (3.90), получим
(Р* - АР* - Р*АТ -Q + Р*£>Р*) Ф2 = 0.
Если 4*2 — невырожденная матрица, то из последнего соотно¬
шения следует, что матрица Р удовлетворяет уравнению (3.89),
т.е. Р* = Р. Но
Ф2 = (Р>Р-ЛТ)Ф2,	Ф2(*0	)	=	Е,
т.е. Ф2 — переходная матрица и, стало быть, det Ф2 Ф 0. Соот¬
ветственно матрица Р может быть найдена по формуле
Р = Ф1Фз1.
Рассмотрим некоторые частные случаи.
1.	Пусть Н = 0. Тогда в (3.89) D = 0, и уравнение Риккати
трансформируется в дисперсионное уравнение. Обозначим че¬
рез Ф(£,£о) переходную матрицу, удовлетворяющую уравнению
Ф(^> ^о) = А(£)Ф(£, to), Ф(*0,*о) = Д.
Легко показать, что
ф2(«) = фтМ,
78

ф1(0	=	$(t,to)^Po	+	J$(to,T)Q(T№T(t0,T)dTj,
P(t) = Ф(Мо) (po + JФ(<о,г)д(г)Фт(*о,г)йт) Фт(г,ь0).
to
2.	Пусть Q = 0. Тогда
Ф1 (t) = Ф(Мо )Po,
Ф 2{t) = Фт(*о,<) + yj ФТ (T,to)D(^(T,to)dTj P0,
P(t) = Ф^1.
3.7.2.	Некоторые условия устойчивости
фильтра Калмана
Обратимся к алгоритму оценки (3.88). Уравнение ошибок оцен¬
ки относительно величины Ах = х — х имеет вид
Ах = Ak Ах + q — Кг, Ак = А — КН.	(3.92)
Важнейшим свидетельством работоспособности алгоритма
(3.88)	является устойчивость (асимптотическая устойчивость)
уравнения ошибок (3.92). Ниже без доказательства приводятся
три теоремы, устанавливающие условия устойчивости фильтра
Калмана для некоторых важных практических случаев.
Теорема. Пусть P(t) — положительно определенная и огра¬
ниченная для всех t > to матрица
(р\Е < P(t) < р\Е, pi = const, р2 = const) ,
являющаяся решением уравнения Риккати
Р = АР + РАТ + Q - PDP.	(3.93)
Тогда уравнение (3.93) устойчиво.
79

Пусть, кроме того, пара (А, Н) наблюдаема. Тогда уравнение
Ах = АкАх	(3.94)
асимптотически устойчиво.
Теорема. Пусть Pi и Р2 — два решения уравнения Риккати
Р = АР + РАТ + Q- PDP, .
соответствующие двум возможным начальным значениям
Pi (to) и P2(t0).
Решения уравнений
Ахг = Ак iAxi,	Aki=A — PiP, (3.95)
Ах2 = Ак2Ах2,	= -А — Р2Р	(3.96)
устойчивы (соответственно асимптотически устойчивы). Тогда
разность АР = Р2 — Р\ ограничена (соответственно стремится
к 0 при t —» 00).
Отметим, что вышеприведенная теорема позволяет положи¬
тельно ответить на вопрос о возможности применения фильтра
Калмана (3.88) в условиях, когда априорные стохастические ха¬
рактеристики динамической системы (3.82) известны приблизи¬
тельно.
Замечание. По условиям теорем требуется ограниченность
матрицы P(t). Установить эту ограниченность можно либо
из физических соображений, либо путем прямого моделирова¬
ния алгоритма оценки.
Рассмотрим стационарный случай: А, Я, Q,P - посто¬
янные матрицы. Поскольку Q >0, существует разложение
Q = ВВТ. Обычно в качестве матрицы В выбирается либо
верхняя, либо нижняя треугольная матрица.
Нас будет интересовать поведение оцениваемой системы при
бесконечном времени наблюдения, а именно условия, при кото¬
рых существует установившееся значение Яоо ? & также свойства
уравнений ошибок
Ах = (А - РооЯ) Ах - Poor + q.
Прежде всего заметим, что при Q = 0 существует решение
Роо = 0 и значит Kqo = 0.
80

Чтобы избежать такого вырождения мы должны, очевидно,
наложить определенные условия на матрицу Q.
Пусть при Q = ВВТ, В такая матрица, что ранг матрицы
W образованной с ее помощью
W= (В,АВ,... ,Ап~1В)
равен размерности системы п.
Матрица W называется матрицей управляемости, а при вы¬
полнении указанного выше условия пара {А, В} называется
управляемой. Понятие управляемости в данной книге не рас¬
сматривается. По этому поводу см., например, [4].
Теорема. Пусть пара (А, Н) наблюдаема, пара (А, В)
управляема. Тогда:
1) Существует единственное установившееся решение Р^
уравнения
Р = АР + РАТ + Q- PDP, Poo = lim P(t),
t—УОО
удовлетворяющее алгебраическому уравнению
АРоо + Р00АТ +Q-РооРРоо = о.
2) Однородное уравнение
Ах = (А- КооН)Ах,	Кос = РооЯтЯ-1
асимптотически устойчиво.
Замечание. В том случае, когда матрица существует и
найдена, уравнение Риккати может быть сведено к линейно¬
му следующей заменой: Р = Р+ + Роо •
Получим для Р*:
Р* = Л.Р» + PtAj - PtDP*,
где А* = А — PooD.
Но последнее уравнение, как мы уже видели, сводится к
линейному путем перехода к обратной матрице.
81

3.8.	Задача сглаживания
(субоптимальный вариант)
Фильтр Калмана применяется для решения задач оценивания в
реальном времени. Это означает, что определение оценок х фа¬
зового вектора х исследуемой линейной системы происходит в
темпе поступления измерительной информации и определяемая
этим алгоритмом оценка x{U) учитывает измерения, получен¬
ные к текущему (реальному) моменту времени U.
Однако, существуют многочисленные приложения, допуска¬
ющие постобработку измерительной информации, что позволя¬
ет иначе сформулировать задачу оценивания. Суть новой по¬
становки состоит в том, что для определения оценки вектора
состояния в каждый момент времени учитывается весь массив
измерительной информации, тем самым используются в соот¬
ветствующем алгоритме как „прошлые“, так и „будущиеа изме¬
рения. Такие задачи носят название задач сглаживания.
Рассматривается стандартная задача оценивания вектора
состояния Xj линейной динамической системы
Xj~|_i = i&jXj -j- Qj,	(3.97)
при помощи измерительной информации zq, z\, ..., z^\
Zj = HjXj + г^,	(3.98)
поступающей в дискретные моменты времени to, t\, ..., tjsf
(tj G
Принимаются стандартные гипотезы для параметров Xq ,
Pq , qj, rj задачи калмановской фильтрации. Требуется в каж¬
дый момент времени tj определить несмещенную линейную
оценку Xj вектора состояния Xj системы (3.97), учитывающую
весь массив измерений zo,	..., zn на интервале решения
задачи [tojtw]- В качестве критерия оптимальности рассматри¬
вается критерий ортогональности. Такая задача носит название
задачи сглаживания на фиксированном интервале.
По-видимому, наиболее конструктивным для численной ре¬
ализации алгоритмов сглаживания на фиксированном интер¬
вале является подход, согласно которому в процедуре опреде¬
ления оптимальной сглаженной оценки xN используются два
82

независимых алгоритма: стандартный фильтр Калмана (пря¬
мой фильтр) и так называемый обратный фильтр - фильтр в
обратном времени [5].
Прямой фильтр Калмана использует „ прошлые “ и „настоя-
щиеи измерения {^о, ..., Zj-i, Zj} и определяет текущую оцен¬
ку:
xf =М [Xj/zo,
При помощи обратного фильтра, использующего „ будущие “
измерения {zj+1, ..., 2/v}, определяется оценка
Хц =M[xj/zj+i,...,zN\.
Искомая оптимальная сглаженная оценка хУ определяется
в результате „склейки“ полученных оценок Xj~ и х^.
Опишем эти алгоритмы.
Прямой фильтр Калмана. Алгоритмы численной реализа¬
ции фильтра Калмана приведены выше и здесь опускаются.
Обратный фильтр [5]. В процессе вычислений, осуществля¬
емых обратным фильтром, определяется оценка щ- вспомога¬
тельной величины, связанной с оценкой х^ исходного вектора
состояния соотношением
Vbj = Ръ^ъу
Здесь Pbj - ковариационная матрица ошибки оценки обратного
фильтра.
Опыт непосредственного применения алгоритма оптималь¬
ного сглаживания показывает, что точная численная реализа¬
ция обратного фильтра и процедуры „склейки“ оценок явля¬
ется плохо обусловленной и может приводить к неустойчивым
вычислениям.
В приложениях хорошо зарекомендовал себя подход, когда
вместо точного оптимального алгоритма, используется субоп-
тимальный алгоритм сглаживания, псубоптимальностьм кото¬
рого, по сути, заключается только в иной инициализации кова¬
риационной матрицы P^N обратного фильтра. А именно, вместо
(PbN) 1 =
83

используется
PbN = ае ■ £, где ае > 1.
Задача оценивания (для обратного фильтра) записывается в об¬
ратном времени, при этом алгоритм ее решения имеет обычную
калмановскую структуру.
Опишем численную реализацию этого алгоритма.
1. Инициализация граничной точки обратного фильтра:
%bN = О’
PbN = ае • Е, (ае - большое число).
2. Этап коррекции:
= %bj + P-bj “ Pj43) »
Pbj = {E — KbjHj)P^,
Kbj = Pb-Hf (HjP^Hj + R3)~l.
Эти соотношения в точности совпадают с соотношениями этапа
коррекции фильтра Калмана в прямом времени.
3. Этап прогноза:
Pbi-x = *7-1 {рьз + Qj-i) (ф7-1)Т •	(3-99)
Соотношения (3.99) структурно совпадают с уравнениями этапа
прогноза прямого фильтра Калмана (3.74). Выделим отличие:
в (3.99) фигурирует обратная переходная матрица Ф”1, в от¬
личие от прямой переходной матрицы Ф в обычном фильтре
Калмана. Матрица Ф-1 может быть определена, при наличии
соответствующей непрерывной модели линейной системы
х = Ах,
например, следующим образом:

4.	„Склейка “ оценок.
Склейку можно осуществить, используя два равноценных
подхода.
1.	Рассмотрим вспомогательную задачу: требуется постро¬
ить оценку х вектора х и определить ковариацию Рах
ошибки оценки Ах = х — х, если заданы два векторных
измерения:
Zj = х + г/, zu = х + Г//,
где
M[ri\ = М[гц] =0, М [r/rf7] = 0,
М [rirJ] =Ri, М [rurfj] = Rn-
Оценка строится как линейная, несмещенная и удовлетво¬
ряющая критерию ортогональности:
х = KjZj + KjjZjj.
Условия несмещенности
К! + Кп = Е,
при этом
Ах = -Kin - Кцгц, zj= Ах + г/, zn= Ах + гп.
Условия ортогональности:
о Т
о г'
м
Axzj
= 0,
м
Axzn
= 0,
отсюда следует
Pax = KjRi, Pax = KuRu
или
Kj = PaxPj j
Кц — РахРц •
85

Из условия несмещенности получаем:
Pax = (RJ1 + Rjiy1 ■
В нашем случае
Zi = x+ , Rj = pf, zji=x^j, Rii = P6“.
Окончательный результат
xf = KjX+ + Kbjx~,
Kj =	Кы	=	р^(рь-у\
pf =
(р/)-1 + (рь-)"1]	(3.100)
2.	Рассмотрим вспомогательную задачу. Оценивается вектор
х по измерению z:
z = х -hr, M[r\ = 0, М [rrT] = R.
Известна априорная информация о векторе х:
fJ>x =	Рх	~~	(*^ А^аг) J •
На основании ранее полученного имеем соотношение
х = fix + K(z - /хх),
где
Рдх = {Е- К)РХ, К = РХ (Рх + R)-1.
В нашем случае
цх = xt, Рх = Р+, г = х“,
Р = Рь“, Рдх=Р/, i =
В результате получаем
if = *++*,(*;,-*;), ifj = ^(p++p-)“I1
86

Р? =	(E-Kj)P+.
Полученные соотношения эквивалентны полученным при
первом подходе, что проверяется непосредственно. Напри¬
мер, величина {Pj*)	,	полученная при первом подходе,
умноженная на Р?* при втором подходе, равна единичной
матрице.
3.9.	Меры оцениваемости
Задача оценивания вектора состояния линейной динамической
системы при помощи доставляемой измерителями информации
приводит к неизбежному вопросу о том, может ли такая задача
быть принципиально решена. Ответ на этот вопрос дает теория
наблюдаемости, позволяющая выделить наблюдаемое подпро¬
странство. Но практически важно знать не только наблюдаема
ли некоторая комбинация компонент вектора состояния, но и
хорошо ли она наблюдаема (или оцениваема). Конечно, если
задача оценивания сформулирована как задача калмановской
фильтрации, то точность оценки полностью определяется ре¬
шением соответствующего ковариационного уравнения.
Однако из указанного решения напрямую не ясно, какую
роль в обеспечении этой точности сыграли измерения, а какую
- априорная информация о системе и ее динамические свойства.
Со многих точек зрения целесообразно строить отдельно ал¬
горитмы для тех переменных состояния, которые могут быть
оценены с удовлетворительной точностью при помощи инфор¬
мации, доставляемой измерителями. Другими словами, нужны
математически формализованные характеристики, определяю¬
щие меры оцениваемости переменных задачи.
С подобной характеристикой мы столкнулись при выясне¬
нии роли собственных чисел грамиана наблюдаемости.
Вопрос о конструктивном введении такой характеристики не
является однозначным. В полном виде она должна учитывать
по крайней мере четыре обстоятельства:
• внутренние свойства системы - свойства пары (Д if);
• диапазон изменения переменных на интервале оценива¬
ния;
87

• время наблюдения (оценивания);
• уровень шумов в измерениях.
Ниже мы рассмотрим некоторые из таких характеристики.
3.9.1.	Сингулярные числа как меры оцениваемости
Обсудим еще одну модификацию алгоритма, доставляющего
оценку по методу наименьших квадратов вектору х в виде ре¬
шения переопределенной системы
z = Нх + r, М [г] = 0, М [ггт] = (т2Е.	(3.101)
Вычислительная процедура, составляющая основу этой мо¬
дификации, в отличие от рекуррентной процедуры, в которой
уточнение оценки происходит от измерения к измерению (см.
п. 3.4.1), имеет дело со всем (полным) массивом измерений с
самого ее начала. Важным достоинством рассматриваемой мо¬
дификации служит получение некоторых положительных чи¬
сел, имеющих смысл мер оцениваемости.
Вернемся к методу наименьших квадратов. Решение пере¬
определенной системы (3.101) в соответствии с этим методом
ищем как решение, минимизирующее функционал J(x):
Известно [12], что любая матрица Н(тхп) ранга к представима
в виде сингулярного разложения
Здесь l/(mxm), V(n х п) - ортогональные матрицы, S(m х п)
- в общем случае прямоугольная матрица и
J(x) = ||Нх — z\\2 = (Нх — z)T (Нх — z).	(3.102)
я = USVT.
f si 0
0 s2 0
0 \
0
S
о
о
О sk 0
о
о
V о
0 /
88

Si > 52 > ... > sk > 0.	(3.103)
Величины Sj (j = 1,2,..., n) называются сингулярными числа¬
ми.
Далее рассматривается частный случай, когда матрица Н
максимального ранга. Тогда
s =	^	^	^	1	Sn	=	diag(s1,s2,...,sn).	(3.104)
Сделаем невырожденные (ортогональные) замены (у - но¬
вый вектор состояния, д - новый вектор измерения):
У = VTX, д = UTz
и коагулируем вектор д:
9 = ( 9п ) ’ ГД6 91^П Х 9п(т ~пх1)-
Вследствие ортогональности использованных преобразований
очевидно, что
min ||Нх — z\|2 = min \\Sy — д\\2.
X	у
Но (см. (3.103))
l|52/-5ll2 = l|5„2/-9/||2 + ||mll2.
Минимум достигается при
У = S-'gu
или, в скалярной форме,
Vj =	7	=	1,2,...,	п.	(3.105)
S3
Найдем ошибку оценки. Имеем
д = Sny + г', г' = £7тг, М [г'г/Т] = сг2Е.
89

Из (3.105) следует:
г'	а2
АУз = Уз ~ Уз = —г> м [АУ2] = Т2 •
bj	bj
Ошибка оценки величины у, порожденная погрешностью ин¬
формации о векторе д, усиливается с коэффициентом усиления
1/sj и, стало быть, число сингулярности может служить коли¬
чественной мерой оцениваемости.
Если сингулярное число мало, то часто оказывается целесо¬
образным в качестве оценки соответствующей переменной при¬
нять нуль и тем самым понизить порядок оцениваемого вектора.
Понижение порядка оцениваемого вектора называется редукци¬
ей.
Приведем без доказательства две важные для приложений
теоремы, во-первых, показывающие, что сингулярные числа
устойчивы к возмущениям матрицы (малые изменения элемен¬
тов матрицы Н приводят к малым изменениям чисел Sj), и,
во-вторых, обосновывающие процедуру редукции.
Теорема. Пусть А, В, Г - матрицы размерности (га х п) и
т > п, Г = А — В. Обозначим сингулярные числа этих матриц,
упорядоченные по невозрастанию, соответственно через aj, /?j,
7j. Тогда
I Pi - Oij\ > Jj.
Теорема. Пусть A - матрица размерности (m x n) и В -
матрица размерности (m х (n — 1)), полученная из А вычер¬
киванием n-го столбца. Если занумеровать сингулярные числа
матрицы В в порядке невозрастания, то они будут разделять
сингулярные числа otj матрицы А следующим образом:
> А > а2 > & > ... > an-i > Pn-i > olu > 0.
3.9.2.	Мера наблюдаемости на основе ортогонализации
Грама-Шмидта
Опишем одну частную меру оцениваемости [2], которая может
использоваться в стационарном случае со скалярным измере¬
нием. Эта мера оказывается полезной для быстрого аналитиче¬
ского выделения мхорошо” наблюдаемых переменных.
90

Итак, рассматривается система
x(t) = Ar(£), z(t) = hTx(t).
Здесь A,h - const и за счет нормировки выполняется условие
\xj(t)\ < 1 на интервале наблюдения.
Введем векторы
д1 =h, д2 = Атд, gj+1 = ATgj.
Матрица наблюдаемости такова
91Т \
9пТ /
Имеют место соотношения
z = xTg\ z = xTg2, z = хт д3,...,	=	хтдп. (3.106)
На базе набора векторов {дг,д2,... ,дп} построим набор орто-
нормированных векторов гх1, ix2, гх3,... , используя процедуру
Г рама-Шмидта:
д1
и1 =	ИЛИ	fl1 = Six«1, «11 = Is1!,
д2 - (д2Т ■ и1) и1
и2 = ——	-f—или д2 = S21U1 + S22U2.
|д2 — (№ • гг) гг|
Точно также получим
Is = ssiu1 + S32U2 + 5331/3.
Если rank N = к, то процесс ортогонализации заканчивается на
к + 1 шаге, когда ик+1 = 0.
Введем новые переменные:

Поведение переменных yi подчиняется уравнениям следующей
структуры:
т = hm + /122/2,
У 2	=	/212/1	+	/222/2 + /23 2/3,
причем в силу соотношений (3.106)
t «22	,	«33
/12	)	/23	, . . . .
5ц	«22
Мерой наблюдаемости \ij переменной 2/j естественно назвать ве¬
личину
= 5Ц/12/2З • • • fj-l,j = 5jj.
Замечание. При практическом использовании введенной здесь
меры нет необходимости прибегать к строгой процедуре ор-
тогонализации. Достаточно только сохранять треугольную
структуру для определения векторов иК Сами эти векторы
следует выбирать так, чтобы углы между ними были доста¬
точно большими, а их модули близкими к единице. Такое до¬
пущение дает возможность при выделении наблюдаемых пе¬
ременных использовать "естественные " комбинации первона¬
чальных переменных, имеющих физический смысл.
3.9.3.	Стохастическая мера оцениваемости
Пусть х = (a?i, #2, ■ • • >£n)T - случайный вектор с априорным
математическим ожиданием рх и ковариацией Рх, {Z} - сово¬
купность измерений, коррелированных с х, х - оценка вектора
х, доставляемая алгоритмом L:
x = L[{Z}}.
Стохастическую меру оцениваемости скалярной величи¬
ны а = сТх, с = (ci, С2,..., сп)Т определим соотношением [7]
 	Я’Да	1	&Аа	л	С?	РдхС
— 1	—	1	\	тП	’	—	№°i	—	^-5
С7ос	аа	у	с1	Рхс
(3.107)
92

где РДх - ковариация ошибки оценки Ах = х — х.
Мера оцениваемости характеризует относительное измене¬
ние среднеквадратичного отклонения сгда ошибки оценки Да
переменной а при ее оценивании посредством алгоритма L.
Для приложений важна задача выделения переменных с ма¬
лыми мерами (если они существуют), которая ставится как за¬
дача отыскания стационарных (экстремальных) значений мер в
зависимости от направления, определяемого вектором с:
(3.108)
Из (3.107) следует, что указанная задача сводится к из¬
вестной задаче нахождения собственных чисел A j собствен¬
ных (главных) векторов регулярного пучка ст (Рдх — АРх) с
квадратичных форм [12]. Пучок называется регулярным, если
Рх > 0.
Характеристическое уравнение
I Рах АРх| = 0
всегда имеет п вещественных корней А(j = 1 , п), которым
соответствуют собственные векторы Ф-7:
PAxV = \jPxVj.
Введем собственную матрицу
ф = (ф1,ф2,...,фп) .
Тогда
ФТРХФ = Е, ФтРДхФ = Л = diag(Ai, А2,..., Ап),
и невырожденное преобразование с = приводит квадратич¬
ные формы стРдхс, стРхс к виду:
Ем? « td
j=1	j=1
Ufl си
дс
= 0.
93

Очевидно, что собственные векторы Ф^ и являются экстремаль-
•т
ными в смысле меры направлениями, а для ctj = Ф-7 х справед¬
ливо
ftoLj ~ 1 x/V
Рассмотрим определение мер оцениваемости в задаче опти¬
мальной калмановской фильтрации.
Непрерывный фильтр Калмана. Матрицы Рх и Рдх нахо¬
дятся из уравнений
Рх	—	АРХ 4- РхАТ “Н Q)
Рах	=	АР&х + РахАт + Q —	РахНтR	^НР&х,
Pxfo)	=	Рдх(*о)«
Дискретный фильтр Калмана. Матрицы Рх и Рдх находят¬
ся из уравнений
Pxj + l = ^jPxj^j + Qj 5
^	=	[E-KjHAP^,
К, = PIxjHj [HjP^HJ + Rj
1 =	+	(3-109)
-1
Для вычисления меры в момент j следует положить PAxj =
Р+
Axj ’
Знание стохастических мер позволяет понять насколько ве¬
лик вклад измерений в оценку тех или иных комбинаций век¬
тора х. В частности, для линейных комбинаций сТх с малыми
мерами может оказаться целесообразной замена их оценок на
априорные значения, т.е. сведение задачи оценивания к мень¬
шей размерности.
Стохастическая мера оцениваемости показывает, в какой степе¬
ни измерения z позволяют улучшить знание компонент вектора
состояния x(t), и является инструментом для выделения комби¬
наций переменных с малыми мерами и возможной последующей
редукции задачи. Но между мерой оцениваемости какой-либо
переменной и ее наблюдаемостью нет прямой зависимости. Из
94

условия «мера равна 1» не следует, вообще говоря, что соот¬
ветствующая переменная наблюдаема.
Пример.
Xi — Х\,	Х2	— Х\ Х2 j
z = х\ + r,	М [r(t)r(s)] = 5(t— s).
Здесь aXl(t) и crX2(t) —> оо при t —> оо, а сгд^^), (Тдхй(^) ограни¬
чены. Отсюда, в частности, следует, что мера /хХ2(£) очевидно
ненаблюдаемой переменной х2 стремится к 1 при t —> оо.
Замечание. При Q = 0 и R = Е уравнение для матрицы, об¬
ратной к ковариационной, Рд* = Я*, в алгоритме калманов-
ской фильтрации то же, что и для грамиана наблюдаемости
Я:
М* = -АТЯ* - Я* А + НТН.
Отличие в начальных условиях
Щ о) = О, M*(to) = P^(to).
Как уже говорилось выше, априорная информация о началь¬
ном состоянии вектора х может интерпретироваться как
дополнительное разовое измерение. В этом случае характери¬
стические числа матрицы Л/**(£) более полно отражают воз¬
можность определения x(t), чем соответствующие характе¬
ристические числа матрицы Я(Ь).
95

4. Калибровка и выставка инерциаль¬
ных навигационных систем
В разделе дается математическое описание двух этапов подго¬
товки ИНС к основному рабочему режиму навигации. Эти эта¬
пы - калибровка (юстировка) и начальная выставка ИНС.
Калибровка чувствительных элементов - необходимый этап
технологического производства навигационных систем, предше¬
ствующий основным режимам эксплуатации ИНС - режимам
начальной выставки и навигации. Калибровка состоит в опреде¬
лении параметров математической модели инструментальных
погрешностей инерциальных датчиков с целью последующей
компенсации этих погрешностей в режиме навигации. Пред¬
полагается, что математическая модель погрешностей задается
априорно.
Начальная выставка ИНС состоит либо в определении на¬
чальной ориентации ее приборного трехгранника относительно
некоторого опорного трехгранника, либо в приведении прибор¬
ного трехгранника в некоторое заданное положение. Выставка
ИНС осуществляется на борту объекта непосредственно перед
режимом навигации.
Решение обеих задач основано на использовании двух при¬
емов: линеаризации и применении алгоритмов калмановской
фильтрации, описанных выше.
4.1. Калибровка ИНС
Под калибровкой понимается определение (оценивание) пара¬
метров инструментальных погрешностей чувствительных эле¬
ментов при помощи проведения экспериментов на специализи¬
рованных стендах.
Традиционные способы калибровки предполагают исполь¬
зование показаний высокоточных датчиков углов стенда и, мо¬
жет быть, датчиков угловых скоростей. При этом обычно нави¬
гационная система последовательно занимает заранее выбран¬
ный набор угловых положений, называемый планом калибров¬
ки. Алгоритмы калибровки строятся с учетом инструменталь¬
ных погрешностей датчиков стенда.
96

При изложении такого метода калибровки требуется доста¬
точно подробное описание конкретного стенда, его характери¬
стик, что в данном пособии не представляется возможным.
Ниже описываются некоторые способы калибровки, либо во¬
все не использующие стенды, либо использующие их таким
образом, что они создают необходимые эволюции приборного
трехгранника ИНС, но при обработке исходной информации по¬
казания датчиков стендов непосредственно не используются [9],
[10].
Замечание. Иногда к калибровке относят задачу определе¬
ния параметров взаимной ориентации приборного трехгран¬
ника ИНС и корпуса объекта, на котором инерциалъная систе¬
ма установлена. В настоящем пособии эта задача, предпола¬
гающая использование дополнительной информации неинерци-
альной природы, не рассматривается.
4.1.1.	Калибровка бескарданной инерциальной
навигационной системы
Бескарданная инерциальная навигационная система (БИНС)
жестко связывается с платформой стенда, которая соверша¬
ет программные вращательные движения в рамках, допускае¬
мых возможностями стенда. Для этого стенд снабжен системой
управления с соответствующими приводами, измерителями уг¬
лов и, может быть, угловых скоростей. Реализация программ¬
ных движений может быть, как это следует из дальнейшего,
достаточно грубой.
В проекциях на оси приборного трехгранника Mz БИНС
измеряются:
• абсолютная угловая скорость и)z приборного трехгранни¬
ка Mz. Результат измерения угловой скорости при помо¬
щи датчиков угловых скоростей (ДУС) обозначается через
• удельная сила fz, действующая на приведенную чувстви¬
тельную массу, расположенную в точке М. Результат из¬
мерения удельной силы при помощи ньютонометров (ак¬
селерометров) обозначается через /' = (/L>/L>/з*)Т-
97

Имеем:
ш'г = uz- uz, /' = fz + Afz.
Знак <—»(минус) у погрешностей ДУС выбран, как это отмеча¬
лось ранее, с целью согласования с обозначениями для систем
с горизонтируемой гироплатформой.
Рассмотрим два варианта решения задачи калибровки
БИНС. В первом варианте основным является понятие при¬
борного трехгранника Mz БИНС, который задается осями чув¬
ствительности ее ньютонометров. Второй вариант основан на
понятии квазиприборного трехгранника Mzx [1], привязанного
к опорному географическому трехграннику.
1.	Используем понятие приборного трехгранника Mz.
Введем модельный трехгранник Му, вектор абсолютной
угловой скорости которого равен сJz. Ориентацию трех¬
гранника Му относительно географического трехгранни¬
ка Мх° определим матрицей ориентации Ly. Имеет место
уравнение:
Ly = ujyLy LyUx.	(4.110)
Здесь их = (0, и cosy?, и sin (р)т - вектор угловой скорости
вращения Земли в осях Мх°, (р - географическая широта
стенда, ujv = oj'z .
Ориентацию приборного трехгранника относительно мо¬
дельного определим вектором малого поворота /Зу
h = {Е + Ру)1у
Справедливо уравнение
Ру = &у@у 4" Vz*
Чтобы замкнуть систему уравнений надо ввести началь¬
ное значение матрицы Ly(to). Вопрос об определении этой
матрицы требует отдельного обсуждения и чтобы не пре¬
рывать связаность текста вернемся к данному вопросу
несколько позже.
98

Введем вектор коррекции
Wy = -g(fz-fy),	/v	=	Ly^0 ).
Имеем
Wy =	^((£?	+	Ду)/у	+	дл-л)	=
— а {Ру^у ^-^г) *
9
Вводя нормированный вектор погрешностей ez:
1 л *
  Д/z?
9
будем иметь
Wy = /3yi3 Н- = —fi/3y + £z,	(4.111)
где Z3 = (Zi3^23) ^зз)Т _ третий столбец матрицы Ly.
2.	Используем понятие приборного квазитрехгранника Mzx.
Квазиприборный трехгранник Mzx определяется соотно¬
шением
lzx = Lylz ,
где I - произвольный вектор.
Ориентацию трехгранника Mzx относительно геогра¬
фического Мх° определим вектором малого поворота
Рх = (/?хЪ Дг2, Асз)т- Легко показать [1], что.
fix — Ly fiy,
а поведение этого вектора подчиняется уравнению
fix = йхРх +	(4.112)
где vx = LyVz.
99

Введем вектор коррекции Wx, равный
WX = I%WV.
Учитывая, что fx = (0,0,#)т, получим
— РхН-
где ех = L^£z.
Или в скалярной форме
Wx\ — ~Рх2 + £х\ч	Wx2	— /3x1 4- £х2> И^З — £хЗ-
(4.113)
Таким образом, в обоих вариантах задача сводится к постро¬
ению оценок вектора состояния линейной динамической систе¬
мы при помощи векторов коррекции Wy или Wx, линейно за¬
висящих от компонент вектора состояния, если математическая
модель инструментальных погрешностей линейно зависит от со¬
вокупности неизвестных параметров, полагаемых константами.
Каждый из таких параметров с* удовлетворяет формирую¬
щему уравнению с* = 0 и вместе с компонентами вектора Ру
или /Зх образует вектор состояния оцениваемой системы.
Как уже неоднократно говорилось ранее, алгоритм оцени¬
вания может быть реализован в двух информационно эквива¬
лентных модификациях: без обратных связей и с их использова¬
нием. Опишем обе модификации, применив для этого частично
операторную форму. Рассмотрим, для определенности, первый
вариант описания задачи калибровки, когда используется при¬
борный трехгранник Mz.
Модификация 1. "Чистое11 оценивание. Для оценок /Зу,
iz имеем соотношения
Ру — &уРу 4- vz -Ь KpAWy,	AWy = Wy +1%РУ — £z,
100

Gxfc, ДНУ = 0, G2[ex,AWy\ = О,
(4.114)
где
Gi^jO] = О, и G20] = О.
Здесь Gi[«] = О и G2\;\ = О - операторы, описывающие форми¬
рующие уравнения, которым удовлетворяют vz и ez.
Коэффициент усиления Кр и операторы Gi, G2 приобрета¬
ют конкретную форму после задания моделей погрешностей и
реализации фильтра Калмана. Обозначим ошибки оценок
Д/Зу = /Зу — {Зу, Avz = Vz — Vzi Абг = Ez ~~ £z-
Тогда уравнение ошибок относительно Д/3 имеет вид
А/Зу = йуА/Зу + Avz - KpAWy,
AWy = Ц A/3y-Aez.	(4.115)
Модификация 2. Оценивание с введением обратных
связей. Введем модельный трехгранник Му*, определив его
матрицей ориентации L* относительно опорного трехгранника
Мх:
L* = u)*L* - L*ux, u*y=ojy- KpAW* - v*z,
AW* = W*- i*z,	W* = !(/' - L*fx),
Gi[v*,AW*] = 0,	G2[e*z,AW*] = 0.	(4.116)
Здесь Кр и Gi, G2 - те же, что и в модификации 1.
Определим рассогласование трехгранников Mz и Му* век¬
тором малого поворота /3*, тогда
/3* = Шур; + Av*z + K^AW*,	(4.117)
где
AW* (3* + Ae*z, Av*z=vz-v*z, Ae*z = ez-s*z.
101

(/3 - третий столбец матрицы L*).
Легко видеть, что
/ = /3*,	^	(4.118)
Для дальнейшего необходимо ввести математическую мо¬
дель инструментальных погрешностей. Будем использовать мо¬
дель, описанную в первой части пособия. Повторим ее описание
и напомним, что она была введена применительно к БИНС с
лазерными ДУС.
Приборный трехгранник Mz введем следующим образом:
ось Mz\ совпадает с осью чувствительности ньютонометра, на¬
зываемого первым; ось Mz<i лежит в плоскости осей чувстви¬
тельности 1-го и 2-го ньютонометров и ортогональна оси Mz\;
ось Mzs дополняет оси Mz\, Mz2 до правого ортогонального
трехгранника.
В соответствии со сказанным, положим
f
{Гп
0
0
е° + Г— + е“,
г = Г21
Г22
0
9
V Гз1
Г* 32
Г33
/ ©и
©12
©13
V*z + Oujy + Vz,
© = ©21
©22
©23
V ©31
©32
©33
(4.119)
Здесь e®, v® - погрешности нулей ньютонометров (нормиро¬
ванные) и ДУС, Tij и 0у, (i,j = 1,2,3i ф j - погрешности
неортогональности осей чувствительности (перекосы), Г** и вц,
(г = 1,2,3) - погрешности масштабных коэффициентов, е®,
- белые шумы с известными интенсивностями.
Компоненты векторов е®, и® и элементы матриц Г, 0 посто¬
янны и подлежат определению.
Как было сказано выше, для каждого из указанных пара¬
метров формирующие уравнения имеют вид
dj = 0.
Введем вектор состояния
С = (А/1, /Зу2, /ЗуЗ, , ^2) ^3 J ©11» ©22, ©33, ©12, ©13, 021, ©23, 031, ©32,
102

£?> £2? £з> Гц, Г22, Гзз, Г21, Г31, Г32) ,
поведение которого определяется ранее приведенными уравне¬
ниями.
Таким образом число оцениваемых величин оказалось рав¬
ным 24.
Величины /3У1, /3У2, Руз могут быть заменены на /?xi, /?Х2, Рхз
с соответствующей заменой уравнений и вектора коррекции.
Об определении начального значения матрицы Ly(to).
Вернемся к вопросу определения начального значения
Ly(to). Возможны различные способы.
Первый способ. Исходим из того, что БИНС устанавливает¬
ся тем или иным способом на платформе стенда, которая грубо
горизонтируется и ориентируется по сторонам света.
Ориентация приборного трехгранника Mz относительно
платформы стенда предполагается также известной с той же
степенью грубости. Оказывается важным только знание на¬
чальной ориентации приборного трехгранника относительно
опорного Мх°.
Под грубостью здесь понимается предположение о началь¬
ном рассогласовании приборного трехгранника Mz с трехгран¬
ником Мх в несколько градусов (например, ~ 5°). Линеариза¬
ция в этом случае еще допустима.
Полагаются заданными априорно начальные значения оце¬
нок вектора состояния с и его матрицы ковариации Рс(^о)-
Именно полагаются: c(io) = 0 и Pc{to) является диагональной
матрицей.
В зависимости от положения приборного трехгранника на
платформе стенда в начальный момент времени матрица Ly(to)
принимает, например, значения
/ 1 0 0 \	/ 0	0	1 \
1. Ly(to) =	0	10,	2.	Ly(to)	=100,
\ 0 0 1 /	\	О	1	о	/
103

О 1 о \
0 0 1.	(4.120)
10 0/
В первом случае модельный трехгранник Му совпадает с
географическим Мх®.
Во втором случае ось Му\ совпадает с осью Мх§, ось Муч
совпадает с осью Мх®, ось Му$ совпадает с осью Мх®.
В третьем - ось Му\ совпадает с осью Мх®, ось Муч совпа¬
дает с осью Мхз, ось Муз совпадает с осью Мх®.
3. Ly(to) —
Второй способ. Приборный трехгранник Mz (платформа
стенда) в течение некоторого времени неподвижен относительно
трехгранника Мх®, жестко связанного с географической коор¬
динатной сеткой. В нашем распоряжении имеется информация
о векторах й и д в виде проекций этих векторов на оси обоих
трехгранников.
Задача может быть решена в пренебрежении инструмен¬
тальными погрешностями чисто алгебраически. На интерва¬
ле "неподвижности" осуществляется осреднение измерений а,
/'. Такое осреднение должно практически устранить высокоча¬
стотные составляющие инструментальных погрешностей, при
этом, разумеется, не устраняются систематические составляю¬
щие.
Далее для осредненных значений оставим прежние обозна¬
чения:
u'z=Wy = (шу1,Ыу2,иуз)Т, /' = (/' 1, f'z2, f'zз)Т.
С точностью до систематических составляющих инструмен¬
тальных погрешностей выполняются соотношения
fz = Lyfx0, Шу = LyUx о	(4.121)
Запишем матрицу Ьу в виде Ly = (1112 13), где
104

Из (4.121) следует, что элементы третьего столбца I3 матрицы
Ьу определяются так:
а также определяется второй столбец I2 матрицы Ly:
7 Myi li3^3
n2	)
U2
где U2 = и cos (р, щ = и sin p.
Первый столбец определим из условия ортогональности ко
второму и третьему, положив его равным векторному произве¬
дению 11 = РТ I3.
Зная структуру матрицы Ly как функцию трех независи¬
мых угловых параметров, например, углов истинного курса ф,
крена 7 и тангажа $, можно при помощи обратных тригономет¬
рических функций определить значения этих углов. Далее все
элементы матрицы Ly перевычисляются, чем обеспечивается ее
ортогональность.
Третий способ. Возможны иные варианты начального опре¬
деления матрицы Ьу при помощи информации /' иш'. Все эти
варианты близки к рассмотренному, но их отличие в том, что
матрица Ly определяется в них сразу как ортогональная. Опи¬
шем один из них.
1. Используем первое из соотношений (см. 4.121). Введем
нормированный вектор и определим столбец 13:
/з = JL
НЛ1Г
2. Используем свойство векторного произведения вектора
угловой скорости вращения Земли на вектор удельной си¬
лы тяжести, взятый с обратным знаком. Результат произ¬
ведения есть вектор, лежащий в горизонтальной плоско¬
сти и направленный на восток.
105

Соответственно введем нормированный вектор	и
II/Z II
/х
И/ХИ
Полученная таким образом матрица Ьу ортогональна, а при
АЛ = 0, vz = 0 удовлетворяет соотношению Ly = Lz.
Найдем ошибки определения матрицы Lz, порожденные по-
грешностями Afz = (Afzl,Afz2,Afz3)T и vz = (yzi,vz2,vz3)T.
Представим матрицу Ly в виде
где AL - линейная часть разложения матрицы Ly в ряд Тейлора
в окрестности величин fz, ujz.
Тогда матрицу ориентации С трехгранника Mz относитель¬
но трехгранника Му можно представить в виде:
С другой стороны — ориентация трехгранника Mz относитель¬
но трехгранника Му определяется вектором малого поворота
/Зу. Имеем
Более наглядное представление о точности определения матри¬
цы L дает вектор (Зх = L%/3y. Имеем
Проводя необходимые выкладки, которые здесь опускаются, по¬
лучим
Ly — Lz + A Z/,
C = E + LZALT.
ЬУ = LZALT.
$x = ALTLZ
или
px = ALTLy.
/3x1
Д/х2
9
106

Ае2
A/xi
АгЗ
(4.122)
где
vx = LyUz, £х = L^£z, £z =
Д/z
9
В рассматриваемом случае, как и ранее, c(to) = 0, а матри¬
ца ковариации Pc(to) не будет теперь диагональной, поскольку
следует учесть соотношения (4.122).
Это предоставляется сделать читателю.
Дополнительная информация для калибровки, связан¬
ная с инвариантностью скалярного произведения двух
векторов к ортогональному преобразованию.
Точность калибровки может быть повышена за счет при¬
влечения некоторых дополнительных измерений. Это особенно
касается систем с достаточно точными датчиками первичной
инерциальной информации. Такой способ получения измерений
основан на идее, изложенной в [а]. Здесь эта идея несколь¬
ко трансформируется. При получении дополнительной инфор¬
мации используется инвариантность скалярного произведения
двух векторов к ортогональному преобразованию.
Рассмотрим временные интервалы, на которых платформа
стенда не вращается, то есть приборный трехгранник Mz непо¬
движен относительно трехгранника Мх°. Предполагается, что
осредненные измерения /', иу (напоминаем, что для осреднен-
ных величин мы оставляем прежние обозначения) содержат
только систематические погрешности.
Сформируем измерения Wg, Wu,	в линейном приближе¬
нии:

(4.123)
Два замечания:
1. В зависимости от ориентации трехгранника Mz относи¬
тельно трехгранника Мх° величины Wg, Wu, W^, очевид¬
но, принимают различные значения.
2. При включении указанных измерений в процедуру кал-
мановского оценивания к ним следует добавлять малые
фиктивные случайные погрешности.
О программных движениях платформы стенда.
Свойства оцениваемости (обусловленности) параметров ка¬
либровки зависят от эволюций установочной платформы, реа¬
лизуемых стендом, или, на информационном языке, от значений
матрицы ориентации Ly(t).
Выбор программного движения платформы стенда - выбор
закона формирования программной угловой скорости Qy(t) -
должен, с одной стороны, осуществляться в рамках его возмож¬
ностей, с другой - обеспечивать достаточно высокую обуслов¬
ленность решения задачи. Такой выбор часто осуществляется
из здравого смысла.
Обозначим через , ППР желаемые (программные) зна¬
чения соответствующих величин. Они определяют программ¬
ный трехгранник Mz11 Р, который в силу неидеальности реа¬
лизации может значительно отличаться от приборного трех¬
гранника Mz. Например, углы рассогласования могут дости¬
гать нескольких градусов.
Рассогласование реального трехгранника Mz с программ¬
ным порождается двумя причинами: несовпадением реальной
угловой скорости платформы П с программной ППР и тем, что
ориентация оси вращения платформы известна, вообще говоря,
не точно. Как это следует из предыдущего, указанное рассогла¬
сование на точность калибровки непосредственно не влияет.
Заметим, что при численном моделировании матрицу Ly(t)
удобнее определять из уравнения
108

О калибровке при помощи одностепенных стендов с го¬
ризонтальной осью вращения.
По числу независимых возможных вращений стенды могут
быть одностепенными, двухстепенными, трехстепенными. Наи¬
более богатые возможности предоставляют последние, в кото¬
рых платформа стенда может совершать три независимых пар¬
циальных вращения.
Но даже на одностепенных стендах можно организовать до¬
статочно большие пространственные эволюции БИНС за счет
ее установки на платформе стенда в разных положениях.
Рассмотрим более подробно калибровку при помощи одно¬
степенных стендов [10]. Существует два основных типа стен¬
дов: первый - ось вращения близка к направлению приборной
вертикали, второй - ось вращения платформы расположена в
плоскости приборного горизонта.
Заранее ясно, что результаты калибровки на стенде перво¬
го типа будут иметь существенно меньшую точность, чем при
калибровке на стенде второго типа. Причина в том, что на стен¬
дах второго типа измеряемая в осях приборного трехгранника
сила / совершает относительно этого трехгранника значитель¬
ные эволюции, в то время как на стендах первого типа вектор /
с точностью до инструментальных погрешностей по отношению
к приборному трехграннику неподвижен.
Вышесказанное подтверждается численным моделировани¬
ем. Далее калибровка БИНС на стендах первого типа не рас¬
сматривается.
Рассмотрим процесс калибровки БИНС для стендов с гори¬
зонтальной осью вращения.
Предварительно обсудим вопрос о выборе (или влиянии)
азимутальной ориентации оси вращения. Два крайних вариан¬
та: первый - ось вращения платформы стенда приблизительно
направлена на восток, второй - ось вращения платформы стен¬
да приблизительно направлена на север.
Умозрительно кажется, что предпочтительным является
второй вариант. Некоторое обоснование: в первом случае при
вращении приборного блока БИНС приведенная восточная со¬
ставляющая систематической погрешности ДУС - восточный
109

дрейф vxi - оказывается постоянной и не наблюдаемой. Уста¬
новлено, что в начальный момент при определении матрицы Lz
наблюдается только комбинация
0x3 ~
и COSip
А при вращении платформы стенда вокруг оси юг-север осу¬
ществляется модуляция величины vx\.
Итак, предполагается, что стенд ориентирован по сторонам
света таким образом, что ось вращения близка к направлению
юг-север и платформа стенда с некоторой степенью точности в
начальный момент горизонтирована. Опишем один из возмож¬
ных алгоритмов калибровки.
Для максимального использования возможностей стенда
(максимально возможной обусловленности решения задачи)
процедура калибровки включает в себя три временных цикла.
БИНС последовательно устанавливается в трех различных по¬
ложениях, определяемых на каждом цикле значениями матри¬
цы Ly.
Для определенности полагаем, что начальное значение мат¬
рицы Ly выбирается по первому способу (4.120), то есть зада¬
ется априорно в соответствии с установкой приборного блока
БИНС на ориентированной по странам света платформе.
Первый цикл: Ly(to) = Е, платформа стенда вращается та¬
ким образом, что
/
cos ф
0
sin^
Ly(t) =
0
1
0
\
. — sin^
0
cos ф
Здесь ijj(t) - угол, на который поворачивается (против часовой
стрелки) платформа стенда.
Второй цикл:
/	0	1 0 \	/0	1 0	\
Ly(to) =	I	0	Oil,	Ly(t) =	j —sinф	0	cosip	I	(4.125)
\	1	0	0/	у cos ф	0	sin ф	)
110

Третий цикл:
О О 1
— sin ф О cos ф
Ly(to) =	1	0	0,	Ly(t)	=	I	cosф 0 sin0 I .(4.126)
Качество калибровки зависит от поведения величины Г2ПР.
Можно показать [10], что наблюдаемость всех переменных с
высокой степенью обусловленности обеспечивается, когда, на¬
пример, каждый цикл включает в себя один или несколько пе¬
риодов, на которых относительная угловая скорость вращения
платформы стенда Г2ПР меняется по закону меандры (прямо¬
угольной "синусоиды”).
Как вариант, на каждом периоде угол ф меняется от 0 до
90° с угловой скоростью П, затем от 90° до —90° с угловой
скоростью — Q и возвращается из положения —90° в нулевое
положение с угловой скоростью П.
Точность калибровки повышается, если на каждом цикле
ввести временной интервал достаточной длительности на кото¬
ром платформа стенда неподвижна. На этом интервале осред-
няются измерения ufz и по результатам осреднения фор¬
мируются величины Wg, Wu, в соответствии с формулами
0 10
0 1 0
(4.123).
В нашем случае на I цикле
Wg  	£z3)
Wu = u2vz2+uzvzZ,
W<p = u2ez2 + usezs — vzz.
(4.127)
На II цикле
Wg = ez2,
Wu = u2vz\ H- usuz2,
111

U2£zl + U3£z2 ~ Vz2.
(4.128)
Ha III цикле
f'z
wa
wu
Wv
Вектор W* = (Wg,WU}Wy,)T в каждом цикле можно рассмат¬
ривать как разовый (одномоментный) вектор коррекции, не со¬
держащий шумовой составляющей. Как уже говорилось, при
его включении в процедуру калмановского оценивания необхо¬
димо введение малых фиктивных шумов. Рассматривая в сово¬
купности величины W* на трех циклах, получаем 9 скалярных
соотношений, связывающих 21 неизвестный параметр.
Начальными оценками и начальными ковариациями для
всех оцениваемых параметров (за исключением составляющих
вектора /3), на каждом цикле служат соответствующие оценки
и ковариации финальные на предыдущем цикле. Для вектора
Р инициализация производится заново независимо от других
циклов. Пример числового моделирования задачи, подтвержда¬
ющего эффективность описанного алгоритма содержится в [10].
Важно подчеркнуть, что рассмотренный способ задания ППР
- только один из многочисленных возможных вариантов, обес¬
печивающих высокую точность калибровки.
4.1.2.	Калибровка ИНС с горизонтируемой
платформой
В описываемом методе используется возможность гироплат¬
формы ИНС поворачиваться в кардановом подвесе относитель¬
но корпуса на достаточно большие углы, так называемые углы
прокачки, причем такие повороты могут быть осуществлены
путем приложения к осям платформы управляющих моментов,
приводящих к вращению платформы с угловой скоростью, про¬
порциональной приложенным моментам. На практике величи¬
= ег1,
= U3VZ1 + U2l/Z3,
=	+	U2£z3	-	Vzl-	(4.129)
112

ны угловых скоростей не должны превышать некоторых задан¬
ных величин. Поворот платформы на достаточно большие углы
позволяет сделать задачу определения параметров инструмен¬
тальных погрешностей хорошо обусловленной (хорошо наблю¬
даемой) .
Методика калибровки состоит в следующем. В начальный
момент приборный трехгранник Mz (трехгранник Mz матери¬
ализован блоком ньютонометров ИНС) горизонтируется и опре¬
деляется его ориентация в азимуте. Для определенности будем
считать, что идеальным положением приборного трехгранника
служит географический трехгранник Мх°, ориентированный
в азимуте в географической координатной сетке. Иная азиму¬
тальная ориентация не меняет сути дела и легко может быть
учтена в алгоритме калибровки. Точность начальной ориента¬
ции приборного трехгранника Mz относительно географическо¬
го Мх° может быть невысокой, лишь бы рассогласование этих
трехгранников могло быть определено вектором малого пово¬
рота.
Рассмотрим более простой (с точки зрения математического
описания) разомкнутый режим калибровки, когда размыкается
канал передачи на датчики моментов управляющего сигнала от
бортового вычислителя ИНС, связанного с режимом навигации.
Командный сигнал управления угловым движением гиро¬
платформы или приборного трехгранника Mz будем формиро¬
вать в виде суммы двух составляющих:
• вектора угловой скорости иу вращения Земли, заданного
в осях модельного трехгранника Му, вводимого ниже;
• управляющего сигнала Пу = (fii, ГЬ) ^з)Т, подлежащего
выбору.
Тогда
Ljz = Q,y Н- иу + vz,
где vz - неконтролируемый гироскопический дрейф платфор¬
мы.
Математическая модель этой погрешности должна быть за¬
дана, ее параметры подлежат определению.
113

Введем модельный трехгранник Му следующим образом.
Его угловая скорость шу равна
со у =	+	иу,
а начальную ориентацию определим единичной матрицей
Ly{to) = Е.
Здесь Ьу - обозначение для матрицы ориентации модельного
трехгранника относительно географического Мх°.
Выбранные начальные условия для матрицы Ly следуют из
предположения, что в начальный момент приборный трехгран¬
ник Mz грубо ориентируется в географической координатной
сетке, то есть горизонтируется и выставляется в азимуте так,
что одна из его осей, а именно ось Mz2, оказывается направлен¬
ной примерно на Север. Другой вариант ориентации - исполь¬
зование режима грубой начальной выставки, описание которой
содержится выше в разделе, посвященном калибровке ВИНС.
Поведение матрицы Ьу подчиняется кинематическому урав¬
нению Пуассона:
Положение приборного трехгранника Mz относительно мо¬
дельного Му определим вектором малого поворота /Зу, так что
lz = (Е + fiy)ly.
Поскольку абсолютные угловые скорости ujz, ujv трехгран¬
ников Mz и Му с точностью до малых второго порядка отли¬
чаются на величину i/z, то уравнение для /Зу примет вид:
Ру = №у + ^у)Ру + vz-
Линеаризованный вектор коррекции Wy построим, исполь¬
зуя
• показания ньютонометров /' = (f'Zl,fZ2,fZ3)T, где
fz = fz + Afz.
Здесь Дfz = (AfZl,AfZ2,AfZa)T - вектор инструменталь¬
ных погрешностей ньютонометров;
114

• Значение fy определяется соотношением
Здесь д - сила, уравновешивающая удельную силу тяже¬
сти, действующая на точку М.
Имеем
Wy = fe-fv) = ^((E + Py)fy + bft-fy) =
= ±(pyfy + Afz).
Вводя обозначения Ly = (1112 Z3), получим
Wy = Pyl3 + ez, где ez =	=	(ebe2, £3)T ■
9
В скалярной форме будем иметь
Wyi	=	А/3 ^23	— Ру2 Ьз + ^1 >
Wy2	=	—ftyzl\3	+ Ад/33 + ^2,
Wys	=	A,2/i3	“ А/1^23 + Сз-	(4.130)
Часто вертикальный	ньютонометр либо не используется,	либо
вовсе отсутствует, и тогда вектор коррекции Wy будет двухком¬
понентным: wy = (wyi, Wy2)T.
Теперь мы можем поставить задачу калибровки более кон¬
кретно. Представим инструментальные погрешности гироплат¬
формы и ньютонометров в виде суммы двух составляющих
vz + vsz и ez+eaz,
где vz, ez - моделируемые составляющие погрешностей; i/J, esz -
немоделируемые составляющие погрешностей, предполагаемые
белыми шумами с известными интенсивностями.
Эти интенсивности обычно определяются заранее при испы¬
таниях по отдельности гироскопов и ньютонометров.
115

Моделируемые составляющие i/z, ez зададим в виде линей¬
ных разложений по известным функциям с неизвестными пара¬
метрами С{. Предполагается, что с* - случайные величины, не
коррелированные между собой с нулевыми средними значени¬
ями и с заданными априорными дисперсиями. Формирующие
уравнения для с* имеют вид: с* = 0. Эти параметры совокупно
с вектором /Зу образуют вектор состояния, поведение которо¬
го подчиняется линейному стохастическому уравнению. Вектор
измерения Wy - линейная комбинация компонент вектора со¬
стояния. В результате оценку величин /Зу, с; можно получить
при помощи стандартной процедуры калмановской фильтра¬
ции.
Для того, чтобы качество оценивания было достаточно вы¬
соким, требуется подобрать управления Пу.
Как пример рассмотрим вариант с двухкомпонентной ИНС,
когда в качестве моделей инструментальных погрешностей ис¬
пользуется модель описанная в [1].
Модель гироскопического дрейфа описывается соотношени¬
ями:
где 0a(i = 1,2,3) - погрешности коэффициентов усиления
в каналах управления датчиками моментов гироплатформы,
Oijihj = 1,2,3) - погрешности ориентации осей прецессии
гироскопов относительно осей приборного трехгранника Мз,
vz —	^2’ uz3)T ~ постоянные составляющие дрейфа vZi ги¬
роплатформы, р = (р1,/02)Рз)Т - вектор отклонения (дебаланс)
приведенной точки подвеса гироплатформы от точки М, задан¬
ный проекциями на оси трехгранника Mz.
Модель погрешности (2-х горизонтальных) ньютонометров
описывается соотношениями:
где £®2 - погрешности нулей, Гц, Г22 - погрешности мас¬
штабных коэффициентов, Г21 - взаимный перекос осей чувстви¬
тельности ньютонометров.
011 012	013
Vz = Ouiy +	+	fzp,	©	=	I	021	022	023
031	032	033
116

В данном случае принимается гипотеза i= 0.
Таким образом, вектор состояния х включает в себя 23 ком¬
поненты:
Pi, 6ц, Oij, Pi,v°i,e°Zk, Tkl, i,j = 1,2,3; к,I = 1,2.
В данной книге мы избегаем расчетов, связанных с число¬
выми характеристиками конкретных приборов.
Заметим, что если рассмотренная процедура калибровки
проводится непосредственно перед включением навигационно¬
го режима ИНС, то она одновременно решает задачу выставки
ИНС (при условии, что конечное положение приборной верти¬
кали близко к направлению географической вертикали). При¬
мер калибровки рассмотренной выше математической модели
инструментальных погрешностей содержится в [9].
В некоторых случаях (например, при больших инструмен¬
тальных погрешностях) целесообразно использовать алгоритм
с обратными связями, так как это было описано в п. 4.1.1.
Приведем этот алгоритм. Но в начале представим для срав¬
нения алгоритм ’’чистого” оценивания, используя для кратко¬
сти изложения в обоих случаях операторную форму. Здесь ис¬
пользуются непрерывная форма фильтра Калмана. Разумеет¬
ся, при численной реализации она должна быть заменена на
дискретную.
Модификация 1. Вариант "чистого" оценивания.
Ру	=	UyPy + vz + KpAWy,
Wy	=	(Wyi,Wy2)T	Wy = (Wyi,Wy2)T,
Wyi = hzPys — 1ззРу2 + £ъ
Wj/2	=	-кзРу3 +кзРУ1 +^2.
A Wyi	=	Wyi-Wyi= l23Ap3-l33Ap2 + Aei + el,
AWy2	=	Wy2-Wy2 = -l13Ap3 + l33Ap1 + As2 + ea2.
Соотношение, доставляющее оценки vz,€z в операторной
форме имеют вид
Gi[P„ AWy] = 0, G2[ez, AWy] = 0,
117

где операторы G\, G2 такие что
GiK,0] = 0, G2[ez, 0] = 0,
описывают поведение величин соответственно vz и ez.
Модификация 2. Вариант оценивания с введением об¬
ратных связей.
Алгоритм обратной связи предполагает введение в сигналы
управления гироплатформой дополнительных составляющих,
зависящих от вектора измерения:
ш*я =шу + V™ + и’ - KfiAW; - й*	(4.131)
Здесь и;* - абсолютная угловая скорость приборного трехгран¬
ника Mz*.
Для абсолютной угловой скорости си* модельного трехгран¬
ника Му* по-прежнему используем
и>1=и>у-	(4.132)
Ориентацию приборного трехгранника Mz*, как результат
управления (4.131) относительно модельного Му*, определим
вектором малого поворота /3*. Уравнение, которому подчиняет¬
ся поведение вектора /3*, имеет вид
$*у = й*ур*у +1/7 + v\ - KpAw; - р;,
Wm =	^*3-/33	P^+eT+el
= li30y3 +	+	e™	+ es2,
aw; = w;-rz.
Величины г/*, e*z определяются операторными соотношениями:
g1[u*z,aw;} = О,	GaK, aw;\ = 0.
В обеих модификациях величина Кр, доставляемая реше¬
нием соответствующих ковариационных уравнений Риккати по
правилам калмановской фильтрации, одна и та же. Информа¬
ционная эквивалентность модификаций следует из легко про¬
веряемых соотношений
АРу = Ру-Ру = 0*,	Vz = v*z, £z = el-
118

Некоторые замечания по поводу построения алгоритмов оце¬
нивания.
1. Обычно применение фильтра Калмана сопровождается
анализом наблюдаемости системы с целью выделения
ненаблюдаемых или слабо наблюдаемых переменных. В
данном случае интерес представляет не то, какие перемен¬
ные являются наблюдаемыми, а то, как наличие оценок
этих переменных (в той или иной комбинации) влияет на
точность решения навигационной задачи. Вместо анали¬
за наблюдаемости проводится анализ поведения элемен¬
тов ковариационной матрицы ошибки оценки в режиме
навигации.
2. Выбор тестовой траектории, на которой проверяется ра¬
бота навигационного алгоритма после того, как в него вве¬
дена компенсация по результатам калибровки, выбирается
из условия, чтобы на ней проявились в наибольшей сте¬
пени все составляющие инструментальных погрешностей.
Выбор управляющего сигнала £ly(t), подаваемого в режи¬
ме калибровки, неоднозначен и определяется тремя обсто¬
ятельствами:
• максимальные углы наклона платформы и угловые
скорости должны удовлетворять заданным техниче¬
ским ограничениям;
• угловые эволюции платформы должны быть макси¬
мально информативны;
• по окончании калибровки платформа должна зани¬
мать исходное положение, за исключением, быть мо¬
жет, азимутального угла.
Точность калибровки может быть повышена за счет при¬
влечения информации, связанной с инвариантностью скаляр¬
ных произведений к ортогональному преобразованию коорди¬
нат так, как это сделано при калибровке Б ИНС.
Этот вопрос здесь не обсуждается, читателю предлагается
рассмотреть его самостоятельно.
119

4.2.	Задача выставки инерциальных
навигационных систем
В соответствии с алгоритмом работы инерциальной навигаци¬
онной системы в бортовой вычислитель должна быть введена
информация о начальной ориентации приборного трехгранни¬
ка. Приборный трехгранник либо ориентируется заданным об¬
разом, либо его положение вычисляется путем обработки дан¬
ных, доставляемых датчиками первичной навигационной ин¬
формации. При этом объект неподвижен относительно Земли
и известны его координаты (например, географические широ¬
та, долгота и высота).
Подобная процедура называется выставкой ИНС на непо¬
движном основании. Здесь рассматриваются только автоном¬
ные методы выставки, т.е. методы с использованием чисто инер¬
циальной информации, и не обсуждается применение таких
средств как теодолиты, астродатчики, магнитные компасы и
т.д.
Замечание. В различных приложениях возникают так назы¬
ваемые задачи выставки ИНС на подвижном основании, ко¬
гда объект, на котором установлена выставляемая система,
движется относительно Земли. При этом источником внеш¬
ней информации для указанной задачи выставки служит ин¬
формация, доставляемая, например, базовым навигационным
комплексом объекта, данные спутниковой навигационной си¬
стемы. В настоящем пособии такие задачи не рассматрива¬
ются.
Автономная выставка основана на том, что датчики первич¬
ной инерциальной информации (ньютонометры и гироскопы)
позволяют определить положение в инерциальном простран¬
стве двух векторов: вектора силы тяжести и вектора угловой
скорости вращения Земли. На базе этих двух векторов может
быть построен тот или иной ортогональный трехгранник. Отсю¬
да, кстати, следует плохая обусловленность задачи автономной
выставки в полярных широтах, когда угол между указанными
векторами мал.
Рассмотрим решение этой задачи применительно к трем си¬
туациям:
120

1. Приборный трехгранник Mz с точностью до инструмен¬
тальных погрешностей неподвижно ориентирован относи¬
тельно инерциального трехгранника 0£.
2. Приборный трехгранник Mz жестко связан с корпусом
объекта (БИНС).
3. Приборный трехгранник Mz жестко связан с горизонти-
руемой гироплатформой.
4.2.1.	Приборный трехгранник неподвижно
ориентирован относительно инерциальной
системы координат
Пусть приборный трехгранник Mz жестко связан с гироско¬
пической платформой в кардановом подвесе и платформа не
управляется. Тогда, если не учитывать дрейфа платформы, она
будет ориентирована неизменно в инерциальном пространстве.
Ориентацию трехгранника Mz относительно трехгранника 0£
определим матрицей Az:
lz —— Azl£ .
Матрица Az подлежит определению. С осями трехгранни¬
ка Mz, с точностью до угловых погрешностей установки, сов¬
падают соответствующие оси чувствительности трех ньютоно¬
метров. Если объект неподвижен, они измеряют составляющие
силы тяжести, взятые с обратным знаком.
Поясним возможность решения задачи выставки. Из-за вра¬
щения Земли вектор силы тяжести описывает в инерциальном
пространстве коническую поверхность, ось которой есть ось
вращения Земли или ось 0£з инерциальной системы координат
0£. Параметры этой поверхности в осях Mz восстанавливаются
по показаниям ньютонометров ИНС. Это позволяет сразу опре¬
делить географическую широту <р места, хотя эта информация
в данной задаче избыточна.
Знание векторов силы тяжести, орта 0£з и долготы Л позво¬
ляет воспроизвести (гринвичский) трехгранник От7, жестко свя¬
занный с Землей. А информация о времени t позволяет опреде¬
лить ориентацию этого трехгранника относительно инерциаль¬
ного трехгранника 0£.
121

Изложим эту задачу на формальном языке. Элементы мат¬
рицы Az являются известными функциями трех независимых
параметров. В качестве таких параметров могут быть выбра¬
ны углы Крылова, Эйлера и т.п. Для определения этих углов,
которые обозначим через	привлечем	уравнения	с	ко¬
эффициентами, зависящими от времени
f'z = -Аг('Ф1,Ф2,Фз)9(-
Здесь fz = fz + A/z измерительная информация, доставляемая
ньютонометрами, причем fz = —gz, a Afz - инструментальная
погрешность ньютонометров.
Справедливы следующие выражения для вектора силы тя¬
жести дг), 9$ в гринвичской и в инерциальной системах коорди¬
нат:
= ~д cos ip cos Л, gv2 = -д cos ip sin Л, дф = -g sin 133)
(cos(Ao + ut)	— sin(Ao + ut)	0	\
sin(A0 + ut)	cos(A0 + ut)	0	I	gv,
о	o	i;
где Ло - начальный угол между осями 0£i	и	Ощ, и - угловая
скорость вращения Земли.
Соотношения (4.133) в каждый фиксированный момент ti
образуют систему из трех уравнений относительно неизвестных
величин ф\, гр2 ,Фз ? но одно из них, очевидно, зависит от двух
других. Рассматривая все их в дискретные моменты с постоян¬
ным шагом At, получим переопределенную систему уравнений,
позволяющую при помощи какого-либо метода осреднения (на¬
пример, метода наименьших квадратов) получить однозначное
решение.
Из-за сложности вычислений (уравнения относительно ^
нелинейны), такой путь нежелателен. Именно поэтому здесь не
приведены явные зависимости элементов матрицы Az от пара¬
метров 'фг. К тому же, чтобы точность определения параметров
фг была достаточно высокой и для устранения влияния дрей¬
фов vz = (^i,i^2, ^з)Т гироплатформы, к соотношениям (4.133)
122

следует добавить математические модели дрейфа vz, и, тем са¬
мым, расширить вектор определяемых параметров.
Предпочтительнее оказывается традиционное двухэтапное
построение решения. Первый этап называется этапом грубой
выставки, второй - этапом точной выставки. Кратко опишем
эти этапы.
Грубая выставка. Грубую выставку можно осуществить, на¬
пример, так. Зафиксируем два момента времени £2- Пока
не будем учитывать инструментальные погрешности в соответ¬
ствующих моделях. Обозначим через /i1^ =	,
/i2) =	\	fzf\ fzf^	векторы составленные из измерений,
доставляемых ньютонометрами в моменты ti, 12. Эти же век¬
торы в проекциях на оси трехгранника обозначим через
4‘> = (Ce/S’f. /f = (С42>.4Т- в«™Рь, /<«,
/(2)
' известны:
= g cos v? cos Ль	= geos ip cos Л2,
f& =	cosy? sin Ab	fj® = geos уз sin A2,
,(i)	.	.(2)
/€з	= 5 sin у?,	= gsiny?,
Л1 = A + Ло + uti, Л2 = A + Ло H-16^2* (4.134)
Далее целесообразно от величин / перейти к нормированным
величинам:
/=-.	(4.135)
9
Введем вспомогательный трехгранник Мст^^з (Мег). Орты
этого трехгранника Oj (j = 1,2,3), заданные своими проекция¬
ми на оси трехгранника Mz, обозначим через aZj, на оси трех¬
гранника - через сг^-.
В инвариантной записи эти орты вводятся следующим об¬
разом:
/(1)	_	/(1)х/(2)	_	гллоа\
<?! = — , 02 = —=	5	,	o3	=	(TiXa2■	(4.136)
\\m и/юх/юц
123

Отсюда следуют соотношения
/*£ \
к =	^2	к,	к	=	k-	(4.137)
В результате получим
lz = ( °z\	&z2	<7z3	)	=	А21$-	(4.138)
To есть матрицу Az можно записать в виде суммы диадных
произведений векторов ajz,
Az = oz i<r|i + az2a\\ + crz3 crj3.	(4.139)
Напомним, что векторы azj определяются по формулам (4.136)
при использовании "идеальных" показаний ньютонометров
/<2), векторы - при помощи этих же формул, но на
основе выражений (4.134) для векторов
Точная выставка. Из-за погрешностей в измерении векторов
g(tj) ньютонометрами и дрейфа приборного трехгранника vz
матрица Az определяется с погрешностью. Вычисленное значе¬
ние матрицы Az обозначим через A!z. Трехгранник с матрицей
ориентации A!z обозначим через Му и назовем его модельным.
Ориентацию трехгранника Mz относительно трехгранника Му
определим вектором малого поворота /3Z:
Задача теперь сводится к построению оценки j}z вектора f3z.
Уравнения измерений ньютонометров, при помощи которых
такая оценка может быть получена, имеют вид
fz = ~9z + A fz = — (^Е + /3Z) 9у + А /г, ду = Аггд$. (4.140)
Перепишем уравнения измерений (4.140) в чуть измененном
виде (компенсируем известную составляющую в показаниях /'
ньютонометров)
^	fz~^~9y	=	-КдУ	4"	Afz — gyPz ~f~ Afz- (4.141)
124

Примем гипотезу постоянства вектора vz. Тогда
Рг = Pz(to)+Vzt,
(4.142)
Примем также гипотезу о том, что погрешности AfZi, (г =
1,2,3) ньютонометров представляют собой независимые белые
шумы заданной интенсивности. Тогда соотношения (4.141), рас¬
сматриваемые в дискретные моменты времени, составляют пе¬
реопределенную систему уравнений относительно шести неиз¬
вестных Pi(to), 02(*о), 03 (*о), »1z, ^2z, V3z-
Ее можно решать, например, методом наименьших квадра¬
тов, либо используя (что предпочтительнее) динамический ал¬
горитм оценивания калмановского типа с применением приемов
регуляризации. В последнем случае вводится вектор состояния
х с компонентами
Формирующие уравнения для его компонент имеют вид
Уравнения измерений (4.141) в скалярной форме таковы
Wi =	01/302	-	03/203	+	A/*l,
W2 = -дузРг + 0з/10з + A/z2>
4.2.2.	Задача выставки БИНС
Задача выставки в рассматриваемом случае состоит в опре¬
делении начальной ориентации приборного трехгранника Mz
БИНС, жестко связанного с блоком ньютонометров и ДУС, а
также в определении некоторых составляющих инструменталь¬
ных погрешностей.
01, 02 , 03, V\z, V2z, V3Z*
Pz = VZ ,	*>2	=	0.
W$ =	0i/20i	-	0з/102	+ A/z3-	(4.143)
Окончательный результат решения задачи - матрица Az
(4.144)
125

Первоначально математически задача ставится так. Прибор¬
ный трехгранник Mz неподвижен относительно вращающейся
Земли и известно точное значение величины tp широты точ¬
ки М. Тремя однокомпонентными ньютонометрами, идеальным
положением осей чувствительности которых служат оси при¬
борного трехгранника, доставляется информация /' о силе fz,
приложенной к точке М:
fz = fz + Afz,
где Afz - вектор инструментальных погрешностей ньютономет¬
ров.
В силу равновесия (точка М неподвижна относительно Зем¬
ли) имеем
fz = ~9zi
где gz - вектор удельной силы тяжести, заданный своими про¬
екциями на оси приборного трехгранника Mz.
Тремя датчиками угловой скорости, оси чувствительности
которых в идеале должны совпадать с осями приборного трех¬
гранника, доставляется информация uz о векторе угловой ско¬
рости ojz, заданном своими проекциями на оси трехгранника
Mz
Vz =шг -vz,
где vz - вектор составляющих инструментальных погрешно¬
стей.
Знак <—» (минус) у погрешностей ДУС выбран, как это
отмечалось ранее, с целью согласования с обозначениями для
систем с горизонтируемой платформой.
Величину ив соответствии с нашими правилами, можно
рассматривать как абсолютную угловую скорость модельного
трехгранника Му (Оу), то есть можно положить u'z = сиу.
Имеем
Uz = Uz,
где uz - вектор угловой скорости Земли, записанный в проек¬
циях на оси трехгранника Mz.
Ранее была введена матрица Lz, определяющая ориентацию
трехгранника Mz относительно трехгранника Мх°, ориентиро¬
ванного в географической координатной сетке
126

Напомним, что их о = (0,гх2,^з)Т? ^2 = и cosip, из = и simp,
Задача выставки БИНС формулируется как задача опреде¬
ления матрицы Lz{tjsi), где - момент перехода от режима
выставки к режиму навигации.
Выставку БИНС целесообразно проводить в два этапа.
Первый этап.
На первом этапе в течении некоторого временного интервала
происходит осреднение величин uj'z, цель которого - умень¬
шить влияние высокочастотных составляющих инструменталь¬
ных погрешностей. Далее матрица Lz определяется алгебраи¬
чески с точностью до систематических инструментальных по¬
грешностей. Соответствующий алгоритм описан в предыдущем
разделе при решении задачи калибровки БИНС (см. п. 4.1.1).
Пусть Uz = Ly - вычисленное значение матрицы Lz. Эта
матрица определяет ориентацию модельного трехгранника Му
относительно трехгранника Мх°.
Ориентация приборного трехгранника Mz относительно мо¬
дельного Му была определена вектором малого поворота /Зу.
Имеет место соотношение:
и компоненты последнего описываются соотношениями (4.122)
/х О = (0,0, у)т.
где I - произвольный вектор.
Было показано, что
Рх = LyPy
/3x1
А/х 2
&х2 j
9
/3x2
Д/х1
9
АеЗ
Vx\ .
	1	£xij
(4.145)
127

где
iyx = LyUZy	£x=Ly£z,
£z = —
9
Таким образом, чисто алгебраически мы определили искомую
матрицу ориентации Lz и ошибки в определении этой матрицы
- вектор Дс, удовлетворяющий соотношениям:
При этом величины ае1,ае2,аез при достаточном времени
предварительного осреднения измерений /' и u)'z, можно поло¬
жить равными нулю (aei = 0, аз2 = 0, аез = 0).
Но на первом этапе выставки также алгебраически можно
определить некоторые параметры инструментальных погреш¬
ностей. Привлечем соотношения (4.123), полученные из усло¬
вия инвариантности скалярного произведения двух векторов к
ортогональному преобразованию систем координат. Они имеют
вид:
— Дс1 + £х2,
= Ае2 — £хЪ
U2 U2
(4.146)
В проекциях на оси трехгранника Мх° получим:
"Wg   £хЗ)
Wu = U2VX 2+U3Vx3,
= U2£x2 + U3£x3 - Vx3.
(4.147)
Откуда следует
ЕхЗ
VX3 ~ U2£x2	=	U3Wg	-	Wp,
128

vX2 + ^з^х2 = — (Wu - u3 ('u3Wg - Wy)).	(4.148)
^2
Обозначим:
ЗВ4	=	ISX2 4” ^3^x2)
*5	—	^x3 — ^2^x2?
ae6	=	£x3-	(4.149)
Окончательный итог первого этапа выставки, который можно
назвать алгебраическим — это определение матрицы Ly и прак¬
тически точное определение величин ае* (г = 1,6) в рамках при¬
нятых гипотез.
Второй этап.
Использование понятия приборного трехгранника
БИНС. Основой математической модели решения задачи
выставки на втором этапе служат уравнения, описывающие
поведение вектора (Зу - вектора малого поворота приборного
трехгранника Mz относительно модельного трехгранника Му.
Модельный трехгранник Му на втором этапе определим как
трехгранник, угловая скорость которого в проекциях на соб¬
ственные оси соу = а;', где lj'z - информация об угловой скоро¬
сти, доставляемая ДУС:
wz=wz- vz.
В качестве модели инструментальных погрешностей vz для
простоты будем использовать, как и на первом этапе, модель,
описываемую соотношением:
Vz=V° + v",
где I/® = 0,	-	белый	шум	с	известной интенсивностью:
му2] = 0, M[vaz(t)vf {s)\ = alE6(t - s).
За начальное положение модельного трехгранника на вто¬
ром этапе принимается то, которое определяется матрицей
129

ориентации Lz(to) - финальный результат первого этапа, т.е.
Уравнения, которым подчиняется поведение (Зу имеют вид:
В случае, когда трехгранник Mz жестко связан с корпусом
объекта, неподвижного относительно Земли, шу заменяется на
иу (сoz — uz). Отметим здесь, что ниже будем рассматривать ва¬
риант определения матрицы Lz, когда приборный трехгранник
может вращаться относительно корпуса объекта с некоторой
заданной угловой скоростью.
Вектор коррекции Wy формируется следующим образом
Ly(t о) — Lz(to).
Поведение матрицы Ly подчиняется уравнению
Ly — (jjyLy Lyuxo.
Принятая модель для величины
Д fz = Af'z-fz = Af°z+Arz.
После нормировки получим:
е° = О, М[еЦ = 0, M[eaz(t)e?(s)} = a2£ES(t-s).
Ру — ШуРу + Vz	1
V°z = 0.
(4.150)
Здесь
(4.151)
130

Представим Ly в коагулированном виде Ly = (I11213). В
результате получим:
Wy = (Зу1 + ez,
где /3 = (/1з,^з^зз)Т-
Задача выставки сведена, таким образом, к построению
оценки вектора Ру при помощи измерений Wy.
Пусть $у оценка вектора /Зу. Тогда оценка Lz находится из
соотношения:
Lz = ^ + P^j Ly.
Использование понятия квазиприборного трехгранника
БИНС. Аналогично тому, как это было сделано при калиб¬
ровке, задачу выставки можно также решать, используя квази-
приборный трехгранник Mzx. Введем векторы /Зж, Wx, vx, ^х,
J.S —	с-0 р9 .
X’
0Х = 1%0у, Wx = LyWy,
vx = vx + v‘ = Lluz, ex = e0x + esx = Llez. (4.152)
В новых переменных получим
Рх = uxoPx + vx.	(4.153)
Запишем соотношения (4.153) в скалярной форме
Рх 1
=
«3/3x2 - «2/3*3 + I'll + ^1,
Рх 2
= '
-«3/3x1 + ^x2 + ^x2.
Рх 3
=
«2/3x1 + VxZ + VxZ,
Wx 1
= — /3x2 + £xl + £xl>
Wx2
= Arl + £°2 + ex2>
wx з
= £x3 + ex3-
Последнее соотношение служит для определения величины £%•
131

Очевидно, что если составляющие векторов v®, постоян¬
ны, величины	также	постоянны.
Из соотношений
M]vsz(t)vf {s)\ =alE-S(t-s), M^eaz(t)ef (s) = cr^E • 5{t-s),
следуют соотношения
M [I/x(*)I/xT(s)] =	■	E5(t	- s), M je*(£)ef (s) = o\ ■ ES(t - s).
Проведем результаты анализа наблюдаемости системы
(4.154) при Lz = const. При анализе используем тот факт, что
если наблюдаема какая-то переменная, то наблюдаема и ее про¬
изводная. Наблюдаемыми оказываются переменные
aei
=
0x1 +
ае2
=
0x2 -
аез
=
о vxl «3 0
Рх з	еХ1,
и2 и2
ае4
=
1/°2 + и3£°х 2,
«б
=
^хЗ ~ и2^х2<
аеб
=
ехЗ-
= 1,6) удовлетворяют соотношениям:
aei
=
u3ae2 - и2ае3 +
ае2
=
-и3 aei + ае4 + i/*2,
«3
=
u2aei + ае5 + i^3,
ЗЭ4
=
0,
=
0,
аев
=
0.
(4.155)
(4.156)
Wx 1	=	—ае2 + е*1,
Wx 2	=	aei	+	е*2,
Wx з	=	зе6 + е*3.	(4.157)
132

Таким образом, задача выставки сведена к построению алго¬
ритма оценивания наблюдаемых переменных ае* при помощи
корректирующего сигнала с составляющими Wx 1, WX2, Wx3.
При принятых моделях инструментальных погрешностей в
качестве такого алгоритма используется фильтр Калмана.
Далее будем рассматривать только вариант обычно имею¬
щий место на практике, когда величинами вида e^Uj можно пре¬
небречь по сравнению со значениями vf. В качестве начальных
условий следует принять соотношения:
ae(io) = 0,	P(t0)	=	[ae(£0)aeT(*o)]	=
/
0
0
0
0
0
0
\
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
<Г1
0
0
0
0
0
0
0
\
0
0
0
0
0
а2е
/
Для устранения вырожденности матрицы Рае(^о) к ней до¬
бавляется регуляризирующая матрица:
P(t0) —> P(to) + А*2
(1
0
0
0
0
0 ^
0
1
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
^0
0
0
0
0
0 )
где /i - малый параметр.
Легко видеть, что при времени оценивания близкой к време¬
ни осреднения измерений на первом этапе, финальные оценки
ае*, (г = 1,6) практически близки к оценкам вектора ае на пер¬
вом этапе (с учетом инвариантных соотношений).
Из сказанного ясно, что алгоритм, доставляющий оценку ае
на втором этапе не улучшает оценку ае, полученную на первом
этапе, когда приборный трехгранник неподвижен относительно
Земли, то есть в рамках используемых гипотез он не нужен.
На практике, однако, во время процесса выставки корпус са¬
молета (или другого объекта), вообще говоря, не является непо¬
движным относительно Земли. Возмущения корпуса могут вы¬
зываться вибрацией из-за работы двигателя, порывами ветра,
133

перемещением внутри самолета экипажа, и тогда применение
указанного алгоритма на II этапе может оказаться оправдан¬
ным. Для конструктивного анализа ситуации с возмущениями
требуется построение более или менее адекватной модели таких
возмущений. Но здесь этот вопрос не рассматривается.
Особо обсудим вопрос об азимутальной ошибке /?хз. Это
i/°
величина содержится в оцениваемой комбинации /Зхз — ■-,
откуда следует, что ошибка определения /Зхз может оказаться
недопустимо большой, особенно в высоких широтах.
Далее величины 1/£1? v^.2,	будем называть соответственно
восточным, северным и азимутальным (вертикальным) приве¬
денными дрейфами ДУС.
Улучшение точности возможно только за счет определения
приведенной восточной составляющей постоянной погрешности
ДУС. Это достигается, если приборный трехгранник будет за¬
нимать разные угловые положения относительно трехгранника
Мх°. В этом случае меняется зависимость величины от ве¬
личин i/^2 > uz3- Для этого должно быть либо предусмотрено
устройство, позволяющее поворачивать приборный трехгран¬
ник относительно корпуса объекта, например, вокруг какой-то
из осей этого трехгранника, либо использовать, например, по¬
ворот корпуса самолета в режиме мрулежки” на аэродроме, при
условии что самолет при этом совершает некоторое число оста¬
новок. Важно только чтобы повороты осуществлялись вокруг
осей, составляющих с осью Мх® достаточно большой угол. Луч¬
ший вариант когда такой угол равен 90°.
Вариант 1. Вращение корпуса БИНС относительно кор¬
пуса объекта.
Корпус объекта неподвижен. Приборный трехгранник БИНС
может поворачиваться относительно этого корпуса. Для опре¬
деленности будем считать, что такой поворот осуществляется
вокруг оси М^з, направление которой близко к направлению
местной вертикали.
Реализуется двухэтапный способ выставки, причем на пер¬
вом этапе задача решается алгебраически, так как это было
рассмотрено выше.
134

На первом этапе определяется матрица Ly(to) = (kj). В
рамках гипотезы о постоянстве погрешностей е®, ошиб¬
ка вычисления матрицы Ly (вектор малого поворота /Зх =
(/3x1, /3x2, /Зхз)Т) такова:
0X1	=	-42,	Рх2	= е°х1, /?хз = -^-,	(4.158)
ucoscp
причем
vl = Ly (*оК,	4 = Ly (to)e°z-	(4.159)
Набор оцениваемых переменных (вектор состояния) можно вы¬
брать различным образом. По-видимому, с учетом возможных
различных ситуаций, удобным является набор, который обозна¬
чим вектором с:
с = (0X1, 0X2, 0X3, V°zl, V% ,^3,6°!, £z2,£z3)T-
Используем соотношения:
0х =	йх0х + LT uz + i/a,
*2 =	о,
= 0.	(4.160)
Вектор коррекции Wx = (Wxi, Wx2, Wx3)T (см.	(4.150), (4.151),
(4.152))
Wx = $х(0,0,1)Т + ЬУ2+е*х =
=	(—0x2,0xi, 0)т + Ly e°z + e%.	(4.161)
В начальный момент второго этапа для ускорения процесса
выставки целесообразно использовать в качестве разовых из¬
мерений инвариантные соотношения Wg^W^W^ (4.146).
Поведение матрицы Ly(t) подчиняется уравнению
Ly(t) = Cjy(t)Ly(t) - Ly(t)uxo(f),	(4.162)
с начальным условием Ly(to), являющимся финальным резуль¬
татом первого этапа выставки.
135

Начальными условиями на втором этапе служат величины:
c(to)	= о,	(4.163)
/ Р00	PPv	РРе	\
Pc(t0) =	Pufi	PuV	Р»е	•	(4.164)
Vxl
«2 '
Р00
РРи
Р0е
Р»0
pvv
Рие
РФ
Реи
Рее
В следствии соотношений
Aei(^o) = £x2i Ап2(^о) = £xl, Асз(^о) =
имеем:
рм =
0
0 \
(0
0
0
0
*2
II
&
0
0
0
2
0
0
4 /
V 0
0
и2
Pp,=
-о2
0
0
0
<те
0
0
0
0
(4.165)
Pv0 — PJV1 Pvv — atP'A, Pep — P0ei Pev — 0, Pee — О^Ез-
Для того, чтобы устранить очевидную вырожденность
Pc(to) следует ввести регуляризирующую добавку, например к
матрице Ррр добавить матрицу /х2Ез, где /х - малый параметр.
Фигурирующая в (4.162) величина ujv представляет собой
вектор абсолютной угловой скорости модельного трехгранника
Му: и)у = иу + Sly, где векторная величина Яу = Clz подле¬
жит выбору. В частности, если вращение осуществляется во¬
круг третьей оси Cl = (0,0, fi)T, и требуется выбрать величину
Cl такой, чтобы обеспечивалась приемлемая точность выстав¬
ки за приемлемое время. Один из возможных вариантов - так
называемое двойное гирокомпасирование.
Циклограмма режима двойного компасирования включает в
себя три цикла и, соответственно, три интервала: (to, Ц), (Ц, £г),
(*2,*з)-
136

Ha I, III циклах (интервалах (to,ti), (£2^3)) корпус БИНС
неподвижен: ujz = uz. На II цикле (интервале (^1,^2)) прибор¬
ный трехгранник БИНС поворачивается относительно корпуса
самолета вокруг оси Mz$ (обычно на угол 90° или 180°) с угло¬
вой скоростью £1. На этом интервале
Wyl = uz 1 — Vz\, UJy2 = uz2 — Vz2, ^3/3 — Uz3 + ^ “ ^z3y
и fi, - величина, подлежащая выбору.
Этот выбор зависит от возможностей поворотного устрой¬
ства.
На всех трех интервалах реализуется описанный выше ал¬
горитм, доставляющий оценку компонентам вектора с.
Отметим, что в рассматриваемой постановке задачи I цикл,
вообще говоря, оказывается лишним. Это означает, что после
этапа "алгебраической” выставки может следовать II цикл ком-
пасирования.
Замечания.
1. На практике выбор начальной ковариации Рдс до опреде¬
ленной степени произволен, то есть допускается нару¬
шение "правильности" задания этой матрицы.
2. При реализации алгоритма оценивания иногда целесооб¬
разно использовать обратные связи подобно тому, как
это неоднократно демонстрировалось ранее.
3. В некоторых случаях для повышения точности оцени¬
вания следует использовать более полную математиче¬
скую модель инструментальных погрешностей. Напри¬
мер, при достаточно большой величине П, возможно,
следует включать в модель погрешности масштабных
коэффициентов ДУС.
В этом случае целесообразно величину Q, задавать в виде
миандры (прямоугольной синусоиды).
4. При малых величинах погрешностей алгоритм оце¬
нивания упрощается за счет пренебрежения этими ве¬
личинами.
137

5.	В алгоритме оценивания на 2-ом этапе вектор /Зх мо¬
жет быть заменен на вектор /Зу. Вектор состояния с в
этом случае имеет вид:
С {fiyl'i Ру2у	^яз)
Вариант 2. Выставка в режиме "рулежка".
Циклограмма выставки в этом случае должна включать в
себя по крайней мере три интервала (£o,£i),	и	(£2,^3)- На
интервале (to, Ц) самолет неподвижен и строится оценка состав¬
ляющих вектора с при помощи стандартного фильтра Калмана.
На интервале (Ц,^) самолет движется по аэродрому. На этом
интервале алгоритм оценивания переходит в режим прогноза.
В нашем случае он описывается уравнениями, которым подчи¬
няются составляющие вектора с
с= (pxl,^x2,0xZ,j>zl>j>z2,t>z3,£°zl^z2^z3) •
Имеем
йхрх +
О,
О,	(4.166)
Ру{1)	—	ujy(t)Ly(t) Ly(t)ux,
РАс = М [ДсДст] .	(4.167)
В непрерывном времени
Рас = АРас + РасАт + Q,	(4.168)
где А - матрица системы (4.166),
q=(^E3x з ОзхбУ	(41б9)
\	^6x3	^бхб	)
Начальные условия Ly(t\), с(Ц), Рдс(Ц) - финальные резуль¬
таты оценивания на интервале (to,ti).
138

На интервале (^Дз) самолет остановлен и оценка вектора
с строится также, как и на интервале (to^i)* Необходимые на¬
чальные условия в этом случае - финальный результат прогно¬
за на интервале
Замечание. На интервалах (£1,^2) и (^гДз) может быть
включен режим навигации. При этом на интервале (^Дз) ис¬
пользуется коррекция по нулевой скорости.
О точном (в рамках принятых гипотез) определение уг¬
ла истинного курса алгебраическим способом.
Как уже говорилось, наиболее трудная часть задачи авто¬
номной выставки БИНС - определение угла истинного курса.
Рассмотренный выше двухэтапный режим выставки обеспечи¬
вает высокую точность и включает в себя как алгебраический,
так и динамические алгоритмы оценивания, использующие кал-
мановскую фильтрацию. Но полученные выше результаты при¬
водят к идее о решении задачи выставки чисто алгебраическим
путем с прежней точностью.
Опишем соответствующий алгоритм. Он включает в себя
так же два этапа.
1. I этап полностью повторяет рассмотренный выше. На этом
этапе после осреднения измерений инерциальных дат¬
чиков вычисляется матрица Ly, определяющая ориента¬
цию модельного трехгранника Му относительно опорного
Мх°. Модельный трехгранник Му - числовой образ при¬
борного трехгранника Mz.
Кроме того, при помощи инвариантных соотношений
(4.147)	(с учетом малости величин e^Uj) оцениваются со¬
ставляющие дрейфа ДУС ихЗ- Здесь	Оцен¬
ки Обозначим Рд;2)^хЗ-
2. На втором этапе приборный трехгранник Mz подворачи¬
вается на некоторый, достаточно большой угол (лучший
вариант - 90°). Поворот осуществляется либо относитель¬
но корпуса самолета при помощи специального поворотно¬
го устройства, либо за счет "рулежки”. Приборный трех-
139

гранник после поворота обозначим через Mz*. Далее по¬
вторяются все процедуры первого этапа:
• осредняются измерения инерциальных датчиков;
• вычисляется матрица ориентации L* модельного
трехгранника Му* относительно опорного Мх°;
• находятся оценки и*2,й*ъ величин	>	гДе	vx	=
Т*Т 77°
Ly VZ•
Имеем
v°x* = L;TLyv°x = Cv°x, с = ь;тьу.
Запишем матрицы Ьу и L* в коагулированном виде:
Тогда для элементов Ctj, i,j = (1,3) матрицы С имеем
сц =
Введем обозначение
V = (йи й2, h)T, V* =	£>2,	йз)т,
здесь	-	оценки величин
Будем иметь: и* = си. Явное решение может быть полу¬
чено несколькими способами, например, с использованием
процедуры МНК.
Аналогично предыдущему, ориентацию трехгранника М*
относительно трехгранника М* определим вектором ма¬
лого поворота /3* и введем вектор малого поворота /3*0:
Тогда L*z = L*(E 4- /3$). С другой стороны аналогично
имеем
R = -еЙ,
&=4и
140

7У°*
/З3* =
U2
В качестве оценок величин Pxh Рх2> РхЗ принимаем вели¬
чины:
Р*1 = -4*2,
Я = 4*1.
7У°*
/О* 	 X1
^3 — —•
U2
В качестве оценки матрицы L* принимаем величину L*,
где
L* = L*{E + j3).
Частный случай.
Поворот осуществляется вокруг оси Mz% на угол, равный
примерно 90°. Матрицу ориентации трехгранника М* относи¬
тельно трехгранника Mz обозначим через Ф. В нашем случае:
/	°	1	°\
Ф =	-1	о	о	.
V	о	о	1 j
Очевидно соотношение:
L* = Ф-L, С — LT^TL,
откуда следует:
^31	~Ьз	+	hlh2	^32	+ ^31?33
С = I h з + I32I31	I32	—hi	+ I32I33
—I32 + h3hi hi + I33I32	^зз
Таким образом для определения величины v\ получаем сле¬
дующую переопределенную систему уравнений:
^2 = (Ьз + h2hl)&l + ^32^2 + (—^31 + 1з21зз)&3\
141

*3 — (—^32 + ^33^31 )^1 + {hi + ^33^32)^2 + ^33^3- (V)
Определим ориентацию трехгранника Mz относительно
Мх° углами курса ф, крена 7 и тангажа v (см. ч. I стр. 123
(5.186)) и предположим, что веичины 7 и v малы. Тогда:
(соБф	йтф	0	\
— Бтф	соъф	О	I	,
о	о	1	/
в этом случае:
V2 = &1,
Vз = ^З-
При этом из уравнения :
р? =	+ (—Z33 + /31^32)^2 + fe + ^31^33)^3
получим: Р* = —V2- Отсюда следует, что в общем случае вели¬
чину v\ целесообразно определять из первого уравнения (V ),
т.е. положить:
^1 = -J	—J—7— (^2 _ ^32^2 + (^31 — ^32^33^3))•
ЬЗЗ + 132*31
Далее очевидным образом определяются величины:
4.2.3. Выставка ИНС с горизонтируемой платформой
В рассматриваемом случае приборный трехгранник Mz жестко
связывается с горизонтируемой управляемой гироплатформой.
Управление осуществляется так, чтобы в идеале он совпадал с
опорным трехгранником Мх, жестко связанным с географиче¬
ской вертикалью и заданным способом ориентированым в ази¬
муте.
Процедура выставки проводится в два этапа. Рассмотрим
алгоритм выставки на первом этапе. Пусть в начальный мо¬
мент платформа тем или иным способом грубо горизонтируется
так, что ориентация плоскости Mz\z<i относительно плоскости
горизонта определяется малыми углами, например, до ~ 10°.
При этом азимутальная ориентация трехгранника Mz произ¬
вольная.
142

Введем вспомогательный трехгранник Мх* такой, что ось
Мхз совпадает с Мхз, то есть задает направление геогра¬
фической вертикали, а ориентация приборного трехгранника
Mz относительно Мх* определяется вектором малого поворота
а = (ai,а2,<*з)Т5 причем аз = 0. Тогда
Информацию об углах ai, а2 доставляют ньютонометры,
оси чувствительности которых с точностью до инструменталь¬
ных погрешностей совпадают с осями Mz\, Mz2.
где A fZl, A fZ2 - инструментальные погрешности.
Чтобы убрать рассогласование а вводится управление (на
датчик моментов по каждой оси подается сигнал, образованный
пропорционально показаниям ньютонометров)
где к - коэффициент усиления, который выбирается допустимо
большим, чтобы удерживать трехгранник Mz в горизонте.
Абсолютная угловая скорость трехгранника Mz в проекци¬
ях на его оси такова:
= Wi + vZl, cuZ2 = W2 + vZ2,	=	W3	+	vzz. (4.171)
где Ws — и sin </?.
Wl = -kff2, W2 = kf[
Отсюда
-kga 1 -kAfZ2 + vZl,
—kga2 + kAfZl +^2,
и sin tp + vZz.
(4.172)
143

После затухания переходных процессов получим
(jjZl = ucos
uZ2 = ucospcosxz,	(4.173)
где Xz ~ азимутальный угол поворота трехгранника Mz вокруг
оси Мхз.
При этом, Xz = Х*> гДе X* угол поворота трехгранника
Мх* относительно трехгранника Мх° вокруг оси Мх3.
Выясним точность решения задачи выставки на первом эта¬
пе. Очевидно, оценки оц, а2 доставляются соотношениями
а,-А.
9	9
Ошибки оценки
Доц = ai - ai = -£z2>	= а2- а2 = £Zl, (4.174)
где ez = (£Zl,£Z2,sZ3)T нормированная погрешность ньютоно¬
метров и
AfZl ^ _ AfZ2
5	сZ2 —
9	9
1
0
-ol2 \
/ COSXz
sinxz
0
0
1
Oil
— sin Xz
COSXz
0
OL2
—Oil
1 /
V 0
0
1
Далее понадобится матрица ориентации Lz трехгранника Mz
относительно Мх°:
Lz =
cos Xz	sm Xz	-a2
— sin Xz	cos Xz	ol	1
ai sin Xz + ol2 cos Xz -<*1 cos Xz + «2 sin Xz 1
Поскольку после затухания переходных процессов трехгран¬
ник Mz неподвижен относительно трехгранника Мх°, то
Ldz — /UZ — LZUX о,
где ихо = (0,и2,из)Т] и2 = ucostp, щ = usmp.
144

Отсюда получим
wZl = и2 sinxz - а2и3,
wZ2 = u2cosx* +<*1U3-	(4.175)
Напомним, что
uZl=W1+vZl> ljZ2 =W2 + vZ2 .
Из (4.175) получим
Щ + i/Zl + a2u3
tgXz =
w2 + vZ2 - aiu3
= ^(1 +
W2 1
Образуем модельный трехгранник Му. Его матрица ориен¬
тации Ly относительно трехгранника Мх° получается из матри¬
цы Lz заменой элементов этой матрицы на модельные значения
Ly
(	cosx7	sinx7	^	^
sin х'	cos х'	^
S-f sinх' - у cosx' -^cosx'- ysinx' 1 j
Числовым образом (модельными значениями) величин cuZl и сoZ2
служат величины Wi, W2. Из этого следует
tsyf=W1+a2ns^mf a2us a1us\
gX W2 - alU3 ~ W2\ Wi W2 )'	(	}
Обозначим через /3Уз = Xz — X* отклонение в азимуте трехгран¬
ника Mz относительно трехгранника Му. Определим эту вели¬
чину. В линейном приближении имеем
tgXz = tg (х' + РУз) ~ tgx' +	•
На основе соотношений (4.176), (4.177), (4.174) и с учетом того,
что с точностью до малых величин
Wi = U2 sin х7,	W2 = и2 cos x7,
145

получим
Руг =	—	(Vzx	cos	х'	-	Vz2	sin	x')	-
U2
- — (e*! cosx' - £z2 sinx') •	(4.178)
Вводя величины
vx\ = yzx cosx “ *^2 sinx , exo = eZl cosx ~ £z2 sinx
получим
A/3
(4.179)
Полученный результат тот же, что и при выставке БИНС.
Как уже говорилось ранее, величина и$ехо для современных
систем на порядок меньше величины vxo. В этом случае вели¬
чину х' можно определить соотношением
щий ориентацию приборного трехгранника Mz относительно
модельного Му в проекциях на оси трехгранника Мх° нахо¬
дится в виде:
Отметим, что в установившемся режиме, после затухания
переходных процессов, трехгранник Mz неподвижен относи¬
тельно трехгранника Мх° и ljz = uz.
Как и в случае БИНС, можно привлечь информацию, ис¬
пользующую инвариантность скалярного произведения двух
векторов к ортогональному преобразованию координат. Сфор¬
мируем измерения Wg, Wu, W^\
В итоге I этапа выставки вектор /Зхо = L%f3y, определяю-
146

Wu =	\	(u2 ~	=
wv> =	~~	{^9 sin <p - fzT^'z)	= \^Vz	~ £*uz-
Проектируя параметры на оси трехгранника Мх°, получим
Wg =	£z3,
Wu = u2vx0 +U3L>xo,
= ux0 - U2exo - u3£x0,
где
vx\ = vZl sin x' + »z2 cos x', vxo = vZ3.
Отсюда находятся величины eZ3, vxo — U2exo, vxo + u$£xo.
Просуммируем итоги первого этапа. Выходными, с учетом
указанных выше пренебрежений, оказываются величины
Ж1	=	Ас° + £х° =	ах\	—	^х? + ех%	— 0)
ж2	Ас§	£х*1	^х§	^х§ ^xj	= 0»
1Ухо
жз =
3 U2
*4	=	^о,
*5 = ^х0=^з-	(4.180)
Перейдем ко второму этапу. Второй этап во многом повто¬
ряет рассмотренный выше этап выставки БИНС. Управление
гироплатформой (приборным трехгранником) будем осуществ¬
лять так, чтобы вектор угловой скорости сoz = (ujZl,bjZ2,(jjZ3)T
приборного трехгранника был следующим
wZl	=	wi	+ vZl,
Uz2	=	CJ2	+ ^z2,
=	W3	+ I/Zs,
где
(jj\ = ucos(psm(x' + Ф),
147

U)2 = W COS cos(x; + Ф),
CJ3 = usin<^-f f2,
и
t
Ф = J П(т)йт,
to
где to - начальный момент второго этапа и одновременно мо¬
мент окончания первого этапа.
Если используются гироплатформы, построенные на базе
механических гироскопов, то высокочастотными составляющи¬
ми v9z можно пренебречь, то есть считать vz = const.
Здесь величина Cl - программная угловая скорость враще¬
ния приборного трехгранника в азимуте, она подлежит выбору.
Ее введение имеет тот же смысл, что и в рассмотренной выше
выставке БИНС.
Выберем в качестве опорного трехгранник Мх*. Матрица
ориентации этого трехгранника Lx+ определяется соотношени¬
ем
/	cos(x' + Ф)	sin(x' +	Ф)	0 \
Lx+ =	[	- sin(x' + Ф)	cos(x' +	Ф)	0 1.
V	о	0	1/
Абсолютная угловая скорость этого трехгранника шх+ =
(Wl,W2,W3)T.
Ориентация приборного трехгранника Mz определяется
вектором малого поворота ах+ = (оц, 0:2, &з)Т- Поведение век¬
тора ах+ подчиняется уравнению
olx* — ^х*(4.181)
По аналогии с тем, как это делалось при выставке БИНС, вве¬
дем вектор ахо = L^ax*. Получим вместо (4.181)
ахо = ихоахо + L%,vz.	(4.182)
Образуем вектор коррекции W* = (PV'f, W£ )Т:
w* = l/' = -a2 + e21+e;1>
148

или в проекциях на оси трехгранника Мх°
(4.183)
где
£1° = е°1 cos(x' + Ф) - e°Z2 sin(x' + Ф),
el° = e°i sin(x' + Ф) + <4 cos(x' + Ф).
К соотношениям (4.182), (4.183) следует добавить уравнения
Вектор состояния с, подлежащий оцениванию, имеет вид
Здесь не приводится анализ наблюдаемости, который по суще¬
ству повторяет анализ наблюдаемости в разделе, посвященном
БИНС. Остаются также в силе все рассуждения этого раздела
относительно определения азимутальной ошибки. Так же как и
в случае БИНС инвариантные соотношения мы будем исполь¬
зовать как разовые измерения в начале второго этапа.
Начальные условия алгоритма оценивания и ковариации на¬
ходятся в соответствии с результатами первого этапа выставки
без учета инвариантных соотношений (см. формулы выставки
БИНС (4.164) и (4.165)).
Один из возможных режимов выставки - двойное гироком-
пасирование. Циклограмма этого режима состоит из трех ин¬
тервалов (t0,ti), (Ц,^), (^2^3)- На интервалах (£о>Ц), (^2^з)
Cl = 0, а на интервале (Ц,£г) Cl = const и чаще всего выбирает¬
ся максимально допустимой. Поворот приборного трехгранни¬
ка на этом интервале вокруг вертикальной оси обычно совер¬
шается на угол 90° или 180°. Кроме того, из-за особенностей
поведения инструментальных погрешностей гироплатформы на
0Zj=0, (j = 1,2,3), ё% =0, (г = 1,2).
С =
(

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

5. Основные виды навигационной
информации, дополнительной
по отношению к инерциальной.
Формирование сигналов коррекции
Для определенности, но не нарушая общности, будем предпо¬
лагать, что в качестве ИНС используется бескарданная инерци-
альная навигационная система (БИНС).
5.1. Информация о высоте
Информация о высоте, как дополнительная по отношению к
инерциальной, уже рассматривалась в первой части пособия
[1]. Было показано, что при помощи такой информации может
быть устранена неустойчивость динамических уравнений оши¬
бок, свойственная автономным инерциальным системам. Рас¬
смотрены некоторые варианты использования для коррекции
БИНС информации о высоте и в частности вариант, в котором
информация о высоте, доставляемая баровысотомером, в виде
обратной связи вводится в "вертикальный канал" ИНС.
Можно показать, что точность навигации по сравнению с
рассмотренными в первой части пособия вариантами повыша¬
ется, если использовать оптимальный алгоритм калмановского
типа.
Сигнал коррекции определяется соотношением:
Wh = ti - ft*, ft* = ft + Р*.	(5.184)
Здесь Ы - модельное значение высоты, ft* - дополнительная
информация о высоте, р* - инструментальная погрешность этой
информации.
Конкретный алгоритм коррекции здесь не обсуждается.
5.2. Угловая информация
В качестве угловой информации рассмотрим сначала информа¬
цию, доставляемую астродатчиками.
151

Положение навигационной звезды в инерциальной системе
<Э£ задается двумя независимыми параметрами, определяющи¬
ми ее орт аё, например, углами возвышения и азимута. В осях
0£ имеем
= (ае«ь0е«2,аеез)Т, IM = 1.
Обозначим через аеу = (aeyi, аеУ2? &уз)т этот орт, записанный в
проекциях на оси модельного трехгранника Оу. Имеем
аеу = А' эе^,
где матрица ориентации модельного трехгранника А' может
быть определена при помощи информации, содержащейся в вы¬
числителе ИНС.
С другой стороны, информация, доставляемая астродатчи¬
ками, позволяет образовывать вектор
ае'г = sez + р*у
где aez = (aezi, эег2> &гз)Т ~ вектор, составленный из проекций
орта аё на оси приборного трехгранника Oz, р* = {pi,р%,Рз)Т
- вектор, составленный из приведенных инструментальных по¬
грешностей астродатчиков.
Здесь мы оставляем в стороне технические детали работы
астросистемы.
Векторы aez, аеу связаны соотношением
aez — (Е -|- /Зу)аеу.
Здесь /Зх - вектор кинематических ошибок ИНС.
Образуем вектор коррекции W*:
W* = ае'г -Э}у = Ру&у + р*у.
Очевидно, что вектор W® определяет (с точностью до ин¬
струментальных погрешностей) только проекцию вектора /Зу на
плоскость, ортогональную вектору аё. Этот факт демонстри¬
руется путем замены величины W^ на информационно экви¬
валентные величины аёуИ^0 и ae^W^0 (векторное и скалярное
произведения орта аеу на вектор коррекции W®).
152

Используя известное правило
а х (Ъ х с) = Ъ(ас) — c(ab),
получим
= Ру-геу(геТ/Зу) + &ур*у,
^yWT = «$Р*у	(5.185)
Одновременное или поочередное визирование двух звезд позво¬
ляет решить задачу определения кинематической ошибки /Зу
полностью. Вопросы построения алгоритмов астрокоррекции
здесь не затрагиваются.
Еще один вид угловой информации: специальными опти¬
ческими приборами измеряется азимутальный угол между на¬
правлением на азимутальный репер (ориентир) и осью Mz2 при¬
борного трехгранника, то есть угол поворота визирной оси оп¬
тического устройства вокруг приборной вертикали.
Обозначим угол между визирной осью и осью Мх2 через ф.
Результат измерения обозначим через ф'. Воспользуемся век¬
торами малых поворотов ах = (ац, аг, Рх = {Р\ч р2> Рз)Т,
7х = (7ъ72>7з)Т> связывающими опорный Мх, квазиприбор-
ный Mzx и квазимодельный Мух трехгранники. Имеем
ф* = ф + Аф + аз.
Здесь Аф - инструментальная погрешность.
Вычислим этот угол при помощи информации ИНС. Резуль¬
тат вычислений
ф* = ф + 7з.
Корректирующий сигнал W^ сформируем в виде
W+ = ф' - ф* = Рз + Аф.
Здесь, как и в случае информации астродатчика, полезным сиг¬
налом измерения является составляющая только кинематиче¬
ской ошибки ИНС.
153

5.3.	Скоростная информация
С точки зрения решения задачи коррекции ИНС по крайней ме¬
ре три датчика доставляют однотипную информацию. Это до-
плеровский измеритель скорости, установленный на самолете;
лаг на морском корабле; одометр - измеритель числа оборотов
колеса, установленный на объекте, движущемся по земле.
Указанные датчики доставляют информацию об относи¬
тельной скорости объекта в проекциях на оси трехгранника,
жестко связанного с корпусом объекта. В свою очередь измеря¬
ются углы между осями этого трехгранника и осями приборного
трехгранника Mz. В случае БИНС вычисляются соответству¬
ющие углы по отношению к квазиприборному трехграннику.
В конечном счете дополнительная скоростная информация
может быть описана векторным соотношением
v; = vz + Pv,
V* = (V1*,Vr2*,Vr3*)T, Vz = (VL,V2,V3)T - вектор относительной
скорости объекта в приборных осях Mz, pv = (pX>pXiP\\)Т ~
вектор приведенных инструментальных погрешностей.
Вектор коррекции Wv описывается соотношением
Wv = V' -V* = SV - ру,
SV - динамическая ошибка определения относительной скоро¬
сти.
В некоторых случаях может отсутствовать информация У3*
о вертикальной скорости.
Еще один вид скоростной информации того же типа добы¬
вается путем периодической остановки движущегося по поверх¬
ности Земли экипажа. На стоянках Vx = 0 и вектор скоростной
коррекции Wv образуется так:
wv = V' = SV.
Таким образом, указанные выше типы скоростной инфор¬
мации не несут напрямую информацию о курсовой (азимуталь¬
ной) ошибке /З3.
154

5.4.	Формирование сигналов корректируемых
ИНС с использованием спутниковой
навигационной информации
Источником позиционной информации до появления спутни¬
ковых навигационных систем (СНС) служили радиосистемы
ближней и дальней навигации (РСБН и РСДН), но в настоящее
время такая информация практически не используется. Ввиду
особой роли спутниковой навигационной информации, которая
содержит в себе как позиционную, так и скоростную информа¬
цию, мы выделяем описание такой информации в отдельный
подраздел.
При построении систем на базе ИНС и СНС обычно исполь¬
зуют термины "комплексные системы" , "интегрированные си¬
стемы" , "тесно связанные системы", "глубоко интегрированные
системы" в зависимости от того, насколько тесно осуществля¬
ется связь этих двух видов навигационной информации. Но со¬
ответствующая терминология в достаточной степени условна.
Ее далее не используем.
Будем различать два вида спутниковой навигационной ин¬
формации: вторичную (выходную) и первичную.
Вторичная информация - это информация о координатах и
скоростях движущегося объекта - в нашем толковании точки М
- в системах координат O77, Ох°, связанных с Землей. Она фор¬
мируется при помощи обработки первичных спутниковых из¬
мерений: кодовых псевдодальностей, доплеровских псевдоско¬
ростей, возможно фазовых измерений. При этом изначально
определяются декартовы координаты г^ и вектор относитель¬
ной скорости объекта в гринвичских осях От], которые за¬
тем пересчитываются в географические координаты Ас, <^с, /гс,
и скорость Vx = (Vf, V2C, V3C)T = (VS, VS, V{})T, где Vg - восточ-
ная, Vjy- - северная, V-- вертикальная составляющие.
Первичная информация - это информация о псевдодально¬
стях и доплеровских псевдоскоростях точки М по отношению
к каждому из визируемых навигационных спутников. Более по¬
дробно об этом ниже.
155

5.4.1. Формирование вектора коррекции при помощи
вторичной позиционной спутниковой
информации
Выходом приемника СНС служат географические координаты
у?с,Ас, hc его антенны. Предположим, что осуществлен соответ¬
ствующий пересчет координат антенны к точке М - приведен¬
ному центру ИНС. Имеем
(РС = (Р + Рф, Ас = А + Рд, h'C = h + p(jl,
где /9д, pch - приведенные погрешности спутниковой инфор¬
мации.
Вектор коррекции	W£)T:
W* =	— (<p' — <pc) = —Acp +
=	(A' — Ac) cos (pc = ДА cos (pc — pcx cos ipc,
Wch = hf-hc = Ah-pch,	(5.186)
где Ap = ^ ДА = A' — A, Ah = Ы — h.
Напомним, что
71*0 = —Aip, 72xo = ДА cos </?c, 73xo = ДА sin </?c,
где вектор 7*0 = (7ix°)72x0)73x°)T - вектор малого поворота
модельного географического трехгранника Му° относительно
трехгранника Мх°.
5.4.2. Формирование вектора коррекции при помощи
вторичной скоростной спутниковой информации
Напомним, что рассматривается ИНС, в которой опорным трех¬
гранником служит трехгранник, жестко связанный с географи¬
ческой вертикалью и заданной ориентацией в азимуте.
При формировании вектора скоростной коррекции следует
учесть, что мы располагаем четырьмя видами информации: по¬
зиционной инерциальной, скоростной инерциальной, позицион¬
ной спутниковой, скоростной спутниковой.
Указанные выше виды информации могут быть представле¬
ны в разных формах:
156

• инерциальная скоростная информация может задаваться
либо проекциями V' = (У/, У2', Уз)т на оси опорного гео¬
графического трехгранника с заданным законом азиму¬
тальной ориентации, либо величинами V#, V^, V{j\
• позиционная и скоростная спутниковая информация мо¬
жет быть отнесена к гринвичской осям Orj:
Vе =	V*	=	(V£,V£,t£)T.
Эта же информация может быть приведена к осям гео¬
графического трехгранника Мх°, определяемого коорди¬
натами (рс, Xе. Результат перепроектировки вектора И? на
оси Мх° - V%,V&y$.
Комбинируя тем или иным способом эти виды информации,
можно получить различные формы вектора коррекции. Рас¬
смотрим три варианта, обычно возникающих на практике.
1.	Скоростная информация СНС задана в гринвич¬
ских осях От/, скоростная информация ИНС в мо¬
дельных осях Оу.
Выходной информацией СНС является вектор	=
относительной скорости в гринвичских
осях От/:
V* = Vv + pcVv,	(5.187)
где ру^ = (pvvi j PvV2 > PvV3)	“ вектор погрешности ско¬
ростной спутниковой информации.
Запишем вектор Vе в проекциях на оси модельного (для
БИНС квазимодельного) трехгранника Оу:
Vyc = ByVvc = Vy + pcVy.
Здесь Ву — В'х - матрица ориентации модельного трех¬
гранника Оу относительно трехгранника От), руу - вектор
погрешности скоростной информации в осях Оу.
Вектор коррекции сформируем в следующем виде:
wi = v;-vyc = Avy-pcVv.
157

Определение величины Д - полной ошибки определения
скорости - дано в I части пособия и
AVv = 6Vy + pyV',
где 5Vy - динамическая скоростная ошибка ИНС, /Зу - ки¬
нематическая.
2. Скоростная информация СНС задана в осях гео¬
графического трехгранника, скоростная информа¬
ция ИНС в модельных осях Оу.
Очевидно, что использование скоростной информации V£o
можно свести к рассмотренному выше случаю, осуществ¬
ляя обратную перепроектировку
v; = (B£ o)Tvxco
и формируя вектор коррекции:
щ = к - уу = К - ВУ (вх0 )т Ксо.
3. Скоростные данные СНС и ИНС заданы в осях
географических трехгранников, ориентированных
в азимуте в географической координатной сетке.
Пусть выходной информацией СНС служат компоненты
вектора относительной скорости V£o = (V£, Vfc, V{j)T. Как
уже отмечалось, вектор V£o представляет собой резуль¬
тат перепроектировки исходного скоростного решения
на оси географического трехгранника Ох°, определяемого
спутниковой позиционной информацией </?с, Ас, hc о гео¬
графических координатах точки М:

Вектор коррекции сформируем в следующем виде:
WZ = V'o-V*o,	(5.189)
где Vyo = (V^, Vjy-, V^) - выходные скоростные данные
ИНС, приведенные к осям модельного географического
трехгранника Оу°:
(cosx'	sinx'	0	\
-sinx'	cosx'	0	1	. (5.190)
0	0	1	У
Х; - модельный азимутальный угол.
Некоторые преобразования составляющих правой части
уравнения (5.189).
(а)
Vyo = S’T(Vy + AVy) = E.,T (vy + 6Vy + pyvf) =
= VyO + Z’T5Vy + Z'TpyS'Z’TV’ =
= VyO + SVyO + fiyO VyO ,
где 6Vyo = (6Ve,SVn,8V3)t, fo = (/3в,Рм,Ра)Т и
SVE
SV" .
Л
@2
cos х;
-sinx'
sinx'
cos х'
cos х'
— sin х;
sin х;
cos х'
in x! )( SVi \
osx' ) \ SV2 )’
(b)
Vxc„ = BcxoB^By0 (vv + pcVn) = Bcx0BT (vyo + pcVy0)
Матрица B^oByO представляет собой матрицу взаим¬
ной ориентации двух числовых образов географиче¬
ского трехгранника Мх°, первый из которых постро¬
ен по спутниковой информации (рс, Ас, второй - по
позиционной информации <рг, А' ИНС.
159

Введем ошибки ИНС по отношению к данным СНС:
Д/ = ц} -	AXd	=	А'	-	Ас.
Как следует из первой части пособия [1], для матри¬
цы В^оВуО справедливо представление:
BcxoB^0=E-jd,
где 7d = (71)72»7з)Т “ вектор малого поворота трех¬
гранника Оу° относительно трехгранника, определя¬
емого матрицей ориентации причем
71 = —A (pd, 72 = AXd cos </?с, 7з = AXd sin ipc.
Тогда
КС° = Vyo-JdVyo+pcVy0.
В результате для измерения Wy (5.189) получим
Щ = Уг -	=	5Vv°	+	(&°	+ 7d) Vy'o - Pvy0- (5.191)
В (5.191) в качестве V*Q можно использовать соответству¬
ющие данные ИНС или СНС.
Особенность представленной модели (5.191) состоит в том,
что величину 7dV^*0 можно вычислить в явном виде по ин¬
формации доставляемой ИНС и СНС, и далее скомпенси¬
ровать. А именно, введем новое измерение Wy :
wf =W£- Yv*0 = SVyO + pyoVy*o - pcVy0.	(5.192)
Перепишем уравнение (5.192) в скалярном виде:
= <5Vi cos x' - 6V2 sin x' + PsVn ~
(/3i sin x' + cos x') Vu ~ Pv0’
V1
Wy2 = 5Vi sin x' + SV2 cos x! ~ РзУе +

(/?! cos x' - z?2 sin x') Vu ~ Pv 0’
v2
wv3 =	6V3	+	(ft sin x' + /З2 cos x'JVg-
(/?i cos x; — #2 sin xO Vtf — py 0. (5.193)
v3
В соотношениях (5.193) в качестве Vg, Vfr, Vfi можно ис¬
пользовать позиционные и скоростные данные ИНС или
СНС.
Модель (5.193) можно упростить, если вертикальная ско¬
рость по крайней мере на порядок меньше горизонталь¬
ной: в первых двух уравнениях следует пренебречь вели¬
чинами, содержащими множитель Уц.
5.4.3.	Формирование вектора коррекции при помощи
первичной позиционной спутниковой
информации
Основой для формирования такого вектора служат величины,
называемые псевдодальностями. Псевдодальность - это резуль¬
тат измерения при помощи СНС расстояния между объектом
- точкой М и j-м навигационным спутником - точкой
Строго говоря, точкой М здесь служит антенна спутниково¬
го приемника. Далее предполагаем, что возможное смещение
приведенного центра БИНС и антенны спутникового приемни¬
ка может быть учтено должным образом.
Введем ряд обозначений:
• с£(М, М^) - расстояние между точками М и М^;
• d(T}, 77W), (1(х,хЩ,	- то же расстояние, выражен¬
ное через координаты точек М и М^\ соответственно, в
осях гринвичской системы координат От/, опорного Ох, и
квазимодельного Оух трехгранников, например:
d (77,77^ = d^’) =	(77 - ,(«)Т (77 — 77W).
Трехгранник Orj жестко связан с вращающейся Землей. Трех¬
гранник Ох жестко связан с географической вертикалью. Ква-
зимодельный трехгранник Оух служит числовым образом трех¬
гранника Ох и определяется информацией, доставляемой ИНС.
161

Имеем
х = Вт],	= Brf*\	у = В'г],	= B'rfj\
где В и В* - соответствующие матрицы ориентации.
Эфемеридная информация, передаваемая с навигационного
спутника, позволяет образовать вектор г]^. Информация о век¬
торе у' и матрице В' содержится в вычислителе ИНС.
Для вектора х и матрицы В = (bij) справедливы представ¬
ления [1]:
х\	=	— е2#яЬ1зЬзз,
%2	=	— е2ДяЬ2зЬзз>
%з	=	Be + h — е2 ИеЬ 33,
/ — sin Л cos х — cos A sin у? sin х — sin <р sin Л sin х + cos A cos х
В = I sin Л sin х — cos Л sin <р cos х — sin Л sin ip cos x — cos Л sin x
у	cos	Л cos <p	sin Л cos ip
cos <p sin x
cos cos x
simp
Из теории CHC известно, что псевдодальность описывается
соотношением (см., например, [6]):
d'(M, M(j)) = Sj)' = d(M, M(j)) + сДт +	,	(5.194)
где с - скорость света; Дт - погрешность часов приемника сиг¬
налов СНС; р^р - приведенная погрешность измерения, в част¬
ности вызванная прохождением радиосигнала через ионосферу
и тропосферу.
Заметим, что термин псевдодальность возник из-за наличия
в выражении (5.194) величины сДт.
Подробное описание принципов работы СНС, алгоритмов,
природы и величин инструментальных погрешностей содер¬
жится во многих источниках и, прежде всего, в публикациях,
посвященным американской системе GPS. В [6], на основе этих
работ, приведено подробное описание математических моделей
и алгоритмов обработки первичных спутниковых измерений
системы GPS.
162

Из (5.194) следует, что определение местоположения точки
М при помощи псевдодальностей требует одновременного визи¬
рования, по крайней мере, четырех навигационных спутников.
Это позволяет путем решения геометрической задачи получить
выходную вторичную (позиционную) информацию и сформи¬
ровать соответствующий вектор коррекции СНС.
Одним из преимуществ коррекции по СНС при помощи пер¬
вичной спутниковой навигационной информации, по сравнению
с коррекцией при помощи вторичной, состоит в том, что прием¬
лемую точность коррекции можно получить при использовании
числа визируемых спутников меньшего, чем четыре.
При формировании вектора коррекции будем использовать
координаты точек М, и модельной точки М' относительно
квазимодельного трехгранника Оух. Корректирующее измере¬
ние имеет вид:
wjp = d'(y, у(й) - d(y',	(5.195)
где
d{y',y(3)) =	=	\J(iу' - у(Я)Т (у' - yU)).
(5.196)
Используя разложение Тейлора относительно модельных пе¬
ременных у' с учетом равенства у' — Ау = у и ограничиваясь
линейными относительно Ау слагаемыми, получим
}	АУ + сАт + Pd]■	(5-197)
В последней формуле величина Ат подлежит оцениванию на¬
ряду с обычным набором ошибок ИНС. В (5.197) приближенно
можно принять у' « (0,0,а)т, где а - длина большой полуоси
навигационного эллипсоида.
Здесь выделим, что модель (5.197) включает ошибку опре¬
деления высоты Sh, поскольку Sh = Дуз. Последнее означа¬
ет, что при построении соответствующих алгоритмов коррекции
используются уравнения ошибок "вертикального канала" ИНС.
163

Опишем известный прием, позволяющий исключить по¬
грешность часов Дт. А именно, один из навигационных спутни¬
ков, например спутник с индексом j (имеющим, как правило,
наибольший угол возвышения), выбирается в качестве базового
и образуются так называемые первые разности измерений:
dli =	,	где	i	=	1,2,..., га, i ф j, .
где m - число видимых спутников.
В этом случае корректирующее измерение не содержит по¬
грешность часов Дт и имеет вид:
wldj = dij' - (d™* - d°>) =
т
где pi = pf - р^\
Замечание. Приведенные выше коррекционные модели позици¬
онных измерений справедливы для случая, когда используются
сигналы измерений одной и той же спутниковой навигацион¬
ной системы ГЛОНАСС или GPS. При совместной обработке
измерений этих систем необходимо дополнительно учитывать
расхождение их шкал времени, что приводит к появлению в
моделях (5.197), (5.198) еще одной неизвестной величины tgps,
характеризующей дробную часть расхождения шкал времени
СНС ГЛОНАСС и GPS. Указанные особенности обсуждаются
в третьей части пособия, посвященного моделям спутниковой
навигации.
5.4.4.	Формирование сигнала коррекции при помощи
первичной скоростной спутниковой информации
Основой для формирования такого сигнала служит радиальная
относительная скорость между точками М и	измеряемая
приемником СНС с использованием эффекта Доплера. Резуль¬
тат измерения обычно называют псевдоскоростью из-за нали¬
чия в измерениях составляющей, зависящей от скорости ухода
часов приемника относительно системной шкалы времени СНС.
164

Как и в случае первичной позиционной спутниковой инфор¬
мации, точкой М здесь служит антенна спутникового приемни¬
ка. Далее предполагаем, что относительное движение приведен¬
ного центра БИНС и антенны спутникового приемника долж¬
ным образом учитывается.
Приемник получает эфемеридную информацию, позволяю¬
щую в реальном времени сформировать координаты и век¬
тор относительной скорости навигационного спутника, от¬
несенных к гринвичской системе координат От).
Обозначим через V= Vd(M, М^) относительную ради¬
альную скорость по линии антенна-навигационный спутник.
Справедливо выражение (для определенности используем век¬
торы, заданные своими проекциями на оси гринвичской систе¬
мы координат Orj):
0)
Vd -	dO)
d^’)	= \J[г] —	(rj — 77^)),	(5.199)
где Vjj - вектор относительной скорости точки М.
Модель доплеровского измерения Vd радиальной скорости
имеет вид:
v}jY = Vjj) + A V^0 + Avjj).	(5.200)
Здесь	-	результат измерения радиальной скорости (5.200)
приемником сигналов СНС; AV^ - общая погрешность допле-
ровских измерений, вызванная уходом часов приемника отно¬
сительно системного спутникового времени; A- приведен¬
ная погрешность измерения, вызванная, в частности, инстру¬
ментальными погрешностями аппаратуры приемника, спутника
и др.
Величину привяжем к модельной системе Оу:
и_ (V-у^)т{у,-Vi1')
d
d (у, yV>)
165

где
yU) = sVJ),	У	=	В'т1,
Kw = B'vy\ Vy = B'V„.
Сформируем оценку Vd^* радиальной скорости	ис-
пользуя выходные (модельные) данные ИНС. Получим'
,а). _
0)
V)-
d0>
Сформируем сигнал коррекции Wy ^ в следующем виде:
w!p = V{dj)' - vjj)*.
Разлагая его в ряд Тейлора относительно модельных пере¬
менных у', Vy с учетом, что
у = у' - Ду, Vy = Vy - AVy,
и оставляя линейную часть разложения, получим:
'.Л т
(j) _ (yij) - у')т
Si)*
=
где
рСЙ* - (Уи)-уу)
SV
Д Vy + pV1 Ay + Д + AV*tj), (5.201)
(t/0) - y') (yO) - y')T
E-
(SV)2
В выражении (5.201) приближенно можно положить
у' « (0,0, а)т.
Так же как и в случае псевдодальностей, можно использо¬
вать первые разности доплеровских измерений:
Kj = (vf' - v?y) - (pf* - Vd0>) ,
позволяющие избавиться от общей погрешности доплеров¬
ских измерений.
166

В результате получим
w(ij) = (.(yw-v')T _ {y(i)-y,)T\ Av +
v	^ dW*	rfO>	)^yv	+
+	(p<*> - P^Y Ay + AV(dij),
AVjij) = AV^-AVP.	(5.202)
Замечание. Модели (5.201), (5.202) содержат члены Р^ТAy,
(р(г) — рО')) Ау, зависящие от позиционной ошибки Ау. Они
могут быть вычислены с использованием текущей оценки Ау,
доставляемой соответствующим алгоритмом обработки. Да¬
лее в измерениях Wy\ Wy^ (5.201), (5.202) разумно осу¬
ществить их алгебраическую компенсацию, что существенно
упростит коррекционные модели первичных скоростных спут¬
никовых измерений.
5.5.	Заключительные замечания
В заключении раздела сделаем несколько замечаний.
После того, как сформирован корректирующий сигнал и вве¬
дена математическая модель инструментальных погрешностей,
построение алгоритма оценивания (фильтра Калмана) стано¬
вится достаточно рутинным делом, хотя при реализации такого
алгоритма в бортовом вычислителе возможны некоторые труд¬
ности.
Первый шаг такого построения - дискретизация исходных
соотношений. При этом алгоритм может быть реализован в
двух вариантах: без введения обратных связей и с обратными
связями.
Построению алгоритма может предшествовать анализ на¬
блюдаемости, позволяющий определить не оцениваемые ком¬
бинации. Но этот анализ более важен для интерпретации ре¬
зультатов математического моделирования задачи коррекции,
стендовых и летных испытаний.
В нестационарных ситуациях анализ наблюдаемости затруд¬
нителен. Но это не столь важно, поскольку большую информа¬
167

цию о качестве алгоритма дают результаты математического
моделирования и результаты летных испытаний.
Более конструктивным является использование стохастиче¬
ской меры оцениваемости, дающее возможность редуцировать
(понижать размерность) задачи. С другой стороны, если в ал¬
горитме оценивания фигурируют ненаблюдаемые переменные,
это не приводит к ошибочным результатам.
Очевидно, что математическое моделирование в реализаци¬
ях не имеет смысла, но оно приобретает смысл, когда привле¬
каются результаты реальных измерений.
Выбор типовых траекторий, на которых проводится моде¬
лирование, осуществляется из условия того, чтобы на этих тра¬
екториях при решении навигационной задачи проявили себя по
возможности все составляющие инструментальных погрешно¬
стей.
168

Литература
1. Голован А.А., Парусников Н.А. Математические основы на¬
вигационных систем. Часть I. Математические модели инер¬
циальной навигации. 3-е издание. М.: Изд-во МАКС Пресс,
2011.
2. Парусников Н.А., Морозов В.М., Борзов В.И. Задача кор¬
рекции в инерциальной навигации. М.: Изд-во МГУ, 1982.
3. Вавилова Н.Б., Голован А.А., Парусников Н.А. К вопро¬
су об информационно эквивалентных схемах в корректи¬
руемых инерциальных навигационных системах. Механика
твердого тела. 2008. N 3.
4. Александров В.В., Болтянский В.Г., Лемак С.С., Парусни¬
ков Н.А., Тихомиров В.М. Оптимальное управление движе¬
нием. М.: Физматлит, 2005.
5. Maybeck P.S. Stochastic Models, Estimation and Control. N.Y.:
Acad. Press, 1979.
6. Вавилова Н.Б., Голован А.А., Парусников H.A., Трубни¬
ков С.А. Математические модели, методы и алгоритмы
обработки измерений спутниковой навигационной системы
GPS. Стандартный режим. 2-е издание. Изд-во Московско¬
го университета, 2009.
7. Варавва В.Г., Голован А.А., Парусников Н.А. О стохастиче¬
ской мере оцениваемости. В сб. "Коррекция в навигацион¬
ных системах и системах ориентации искусственных спут¬
ников Земли". Изд-во МГУ, 1987.
8. Лоусон Ч., Хенсон Р. Численное решение задач метода наи¬
меньших квадратов. М.: Наука, 1986.
9. Парусников Н.А., Тихомиров В.В., Трубников С.А.
Определение инструментальных погрешностей инерци¬
альных навигационных систем на неподвижном основа-
нии//Фундаментальная и прикладная математика. 2005. Т.
11.	№. С. 159-166.
10. Вавилова Н.Б., Парусников Н.А.; Сазонов И.Ю. Калибров¬
ка бескарданных инерциальных навигационных систем при
помощи грубых одностепенных стендов. В сб. Современные
проблемы математики и механики. Том I. Прикладные ис¬
следования. - М.: Изд-во МГУ, 2009. стр. 212-223.
169

11. Измайлов Е.А., Jlene С.Н., Молчанов А.В., Поликов-
ский Е.Ф. Скалярный способ калибровки и балансировки
бесплатформенных инерциальных навигационных систем.
Сборник материалов юбилейной XV Санкт-Петербургской
конференции по интегрированным навигационным систе¬
мам. С-Пб., 2008.
12. Гантмахер Ф.Р. Теория матриц. М.: Наука, 1988.
170

Научное издание
ГОЛОВАН Андрей Андреевич
ПАРУСНИКОВ Николай Алексеевич
МАТЕМАТИЧЕСКИЕ ОСНОВЫ
НАВИГАЦИОННЫХ СИСТЕМ
Часть II
Приложения методов оптимального оценивания
к задачам навигации
2-е издание, исправленное и дополненное
Оригинал-макет изготовлен
издательской группой механико-математического факультета МГУ
Технический редактор Ж.Г. Гаврилова
В оформлении обложки использована фотография Н.Н. Молчанова
Подписано в печать 08.10.2012 г.
Печать офсетная. Бумага офсетная.
Формат 60x90 1/16. Усл.печ.л. 10,75. Тираж 200 экз. Изд. № 385.
Издательство ООО “МАКС Пресс”
Лицензия ИД N 00510 от 01.12.99 г.
119992, ГСП-2, Москва, Ленинские горы,
МГУ им. М.В. Ломоносова, 2-й учебный корпус, 527 к.
Тел. 939-3890, 939-3891. Тел./Факс 939-3891.
Напечатано с готового оригинал-макета
Типография МГУ
119991, ГСП-1, Ленинские горы, д. 1, стр. 15
Заказ 1237.