Разработка и исследование характеристик платформенной инерциальной навигационной системы полуаналитического типа
Введение
Цель настоящей работы разработать алгоритм
платформенной инерциальной навигационной системы, работающей в геоцентрической
системе координат, и определяющей в этой системе следующие параметры:
Координаты
Скорости
Углы ориентации
Так же предусмотрена задача получения
позиционной и скоростной информации в географической системе координат (φ,
λ, 
),
вычисление углов курса, крена и тангажа.
Исследование точностных характеристик системы по
уравнениям ошибок и оценка влияния ошибок начальной выставки и гироскопов на
точность ИНС производится при помощи алгоритма разработанного в программной
среде matlab.
В конструкции ИНС используется акселерометры
А-17 и лазерные гироскопы ГЛ-1 производства Раменского приборостроительного
завода.
Исходные данные
Траекторные условия: полет с постоянной
скоростью W=900км/ч, на постоянной высоте H=10000м, курс постоянный
произвольный, время полета 1.5 часа.
Точностные характеристики системы: дрейфы
гироскопов 0.005-0.05 град/час, начальные ошибки координат 15 м, скорости 0.1
м/с(до 5 м/с), ошибка измерений акселерометра 0,05g.
Краткое изложение теоретических
сведений систем координат, в которой работает представленная ИНС
Геоцентрическая система координат
Рис. 1
Геоцентрическими координатами точки O1 являются:
геоцентрический радиус R, угол 
между плоскостью
экватора и радиусом R, угол λ между
плоскостью, содержащей ось 0ζ и
точку 01 и плоскостью 0ζξ. Пример
геоцентрической системы координат предсатавлен на рисунке 1.
Географическая система координат
Рис. 2
Свяжем с земным эллипсоидом правую ортогональную
систему координат Oξηζ (рисунок
2), при этом начало О совместимо с центром Земли, ось Оζ
напрвим
по малой оси эллипсоида в сторону северного полиса, оси Оξ,
Оη
расположим
в плоскости экватора, причем Оξ - по
линии пересечения гринвичского меридиана с экватором.Возбмем некоторую точку О1
в системе координат Oξηζ
и проведем через нее нормаль к земному эллипсоиду. Положение точки О1 в системе
координат Oξηζ можно
определить углом φ, составляемым
указанной нормалью с плоскостью экватора, углом λ, образуемого
плоскостями меридиана точки О1 и гринвического меридиана, и отрезком отрезком h
от точки пересечения нормали эллипсоида до точки О1. Данные углы φ
и λ
называют
соответственно географической или геодезической широтой и долготой. Величина
отрезка нормали h с большой точностью совпадает с величиной высоты точки О1 над
уровнем океана. Геоцентрическая долгота, очевидно, равна географической.
Алгоритм работы ИНС
Введем систему координат Oxyz с началом в центре
Земли О и с ориентацией одноименных осей по осям платформы (акселерометров).
Выполнение условия, чтобы все время ось Oz платформы совпадала с вектором
положения R, для введенной системы координат Oxyz означает, что
= Rz; x = y = 0,(1)
т.е. ориентация по вектору положения имеет
место, если определяемые координаты х и у равны нулю. Этим и определяется
зависимость ориентации платформы от определяемых координат х, у.
Дифференцируя (1), получим выражение для
скорости:
(2)
откуда имеем:
(3)
Первые два выражения (3) определяют законы
управления ориентацией платформы (измерительных осей акселерометров), т.е.
значения угловых скоростей поворота платформы 
в
функции времени, при идеальной реализации которых выполняется условие (1) и,
таким образом, осуществляется заданная ориентация платформы, т.е. ориентация по
вектору положения R. Что же касается ориентации в азимуте (геоцентрическом
горизонте), то она может выбираться независимо от выполнения условий ориентации
по вектору положения (1).
Если учесть условие (1) и соответственно этому
считать, что осуществляется идеальная реализация законов управления, т.е.
и учесть выражения (3), то получим:
(4)
При различных способах ориентации платформы в
азимуте вид уравнений зависит от закона управления этой ориентацией, т.е. от 
.
Рассмотрим сначала уравнения (4) при ориентации
платформы в азимуте по координатным осям хк, ун сферической системы координат.
Для приведения уравнений (4) к виду, при котором
определяются и две другие сферические координаты Ф и Л, надо знать законы
управления полной ориентацией платформы, т.е. по вектору положения и в азимуте.
Поскольку первые уже известны и соответствуют первым двум выражениям (7.62), то
необходимо установить закон управления ориентацией в азимуте, т.е. ,
реализацию которого должна обеспечить система управления. Затем эти законы
управления ориентацией необходимо связать с производными координат
(5)
Получим выражения проекций абсолютной скорости
вращения координатного трехгранника на оси Ox, Оу, Oz через указанные производные
(6)
Из второго и третьего равенств (5) определится
соотношение
(7)
подставив в которое правую часть первого
равенства (3), получим выражение
(8)
определяющее собой закон управления 
ориентацией
платформы в азимуте в функции времени. При идеальной реализации закона 
,
согласно (6), осуществляется заданная ориентация осей платформы по координатным
осям Охк, Oyк.
Если теперь при идеальной реализации законов
управления ориентацией правые части первого и второго равенств (3) подставить
соответственно в (5), то вместе с (6) получим искомые соотношения
(9)
Используя (7), представим уравнения
функционирования (4) для рассматриваемой ИНС с управляемой ориентацией
трехгранника измерительных осей для случая ориентации по координатным осям
сферической системы координат, считая поле тяготения сферическим и реализацию
законов управления идеальной:
(10)
Представим теперь уравнение (4) для ИНС с
азимутально-свободной ориентацией платформы. В этом случае платформа и
материализуемый ею трехгранник измерительных осей не вращаются вокруг оси Oz по
отношению к инерциальной системе координат. И в данном случае для преобразования
уравнений (4) к виду, при котором определяются также и две другие координаты Ф
и Л, надо знать законы управления полной ориентацией платформы. Закон
управления ориентацией платформы в азимуте в этом случае сводится к ее
стабилизации в азимуте, т.е.
(11)
а законы управления по вектору положения
определяются по-прежнему первыми двумя равенствами (3).
При идеальной реализации законов управления
указанной ориентацией соотношения, определяющие эти законы, необходимо связать
с производными сферических координат Ф, Л. Так как законы управления
реализуются вращением платформы вокруг ее осей, а производные Ф и Л есть
составляющие угловой скорости вращения трехгранника координатных осей, то для
установления связи этих скоростей надо спроектировать компоненты линейной
скорости, определяемых системой, на координатные оси. Пусть оси платформы Ох и
Оу составляют с соответствующими координатными осями Охк и Оук угол 
(рис.
), тогда получим
(12)
Если (7) абсолютная угловая скорость вращения
координатного трехгранника 
вокруг оси Оzк, то
угловая скорость вращения платформы вокруг Оzк по отношению к указанному
трехграннику выразится
(13)
В свою очередь величина 
определится
первым равенством (9), что в новом обозначении запишется
(14)
На основании равенства (11), второго равенства
(9) и (13), (14) получим
(15)
В соответствии с (11), (15) уравнения
функционирования (4) для случая ИНС с азимутально-свободной ориентацией
платформы, считая поле тяготения сферическим, а реализацию законов управления
идеальной, будут иметь вид
(16)
Для перехода к относительным значениям Uв и Uc
необходимо ввести параметры от переносного движения (вращения Земли), после
чего получим
(17)
Ориентация платформы с установленными на ней
инерциальными элементами (акселерометрами, лазерными гироскопами) ИНС
сферической системы координат, для которой выведены уравнения функционирования,
реализуется при помощи управляемых силовых или индикаторно-силовых
гиростабилизированных платформ, а также при помощи управляемой платформы
вращением по отношению к свободной стабилизированной платформе.
Пересчет координат из
геоцентрической в географическую систему координат
Выразим геоцентрический радиус R точки О1 через
модуль вектора земного эллипсоида R1, отрезок h' продолжения этого вектора до
точки O1 и широту ф1. Используя уравнение эллипса в полярных координатах, т.е.
И выражение квадрата экстцентриситета, получим
Получим связь координат 
и
R c 
и
h. Согласно рисунку 1 выразим координаты 
через
R, ,
:
ζ
Получим равенство
На основании которого, используя (1.6), получим
искомое соотношение
Зависимость R от получится,
если взять равенство
И подставить в него выражения для координат
согласно (1.6):
Образуем на сфере с геоцентрическим радиусом R
сопровождающий трехгранник 01х2y2z2, связанный с точкой 01 подобно тому, как
был введен сопровождающий трехгранник 01х1y1z1 поверхности h = const. Ось O1z2
направим по геоцентрическому вектору, ось O1y2 расположим в плоскости меридиана
точки 01 и направим в сторону северного полюса, ось O1x2 направляется так, что
образуется правый ортогональный трехгранник. Ориентация трехгранника 01х2y2z2
по отношению к системе 
определяется
таблицей направляющих косинусов.
Из сравнения трехгранников 01х2y2z2 и 01х1y1z1
видно, что их оси 01х2 и 01х1 совпадают. Данные трехгранники повернуты вокруг
совпадающих осей относительно друг друга на угол ф-ф1, т.е. на величину
разности географической и геоцентрической широты. Взаимное расположение
трехгранников определяется таблицей направляющих косинусов:
y2 z20 0;cos (ф ф1) -sin (ф ф1);(1.20)sin (ф ф1)
cos (ф ф1).
Выражение для разности (ф ф1) определится:
Вследствие малости величин 
и
,
считая также величину h/a малой и раскладывая правые части указанных формул в
ряды по степеням 
и h/a, будем иметь
Ввиду малости и
упрощается
матрица направляющих косинусов. Принимая cos ()
= 1, sin 
)
= 
получим
y2 z20 0;1 - 
1
При подстановке значения =
0,0067 получаем максимальное отклонение истинной вертикали от геоцентрической,
равное 
=
0,00335, что соответствует 
и имеет место на
широте 
.
С увеличением h эта разность убывает, но убывание происходит медленно.
Например, при h = 100 км разность составляет 
.
По этой причине при небольших значениях h можно считать
Анализ ошибок ИНС
Чтобы оценить точностные характеристики системы
воспользуемся точным уравнением, описывающее ошибки ИНС:
(1)
где 
После ряда преобразований и упрощений уравнение
ошибок ИНС, в которой используется внешняя информация о высоте полета, примет
вид системы, состоящей из двух векторных уравнений:
(2)
Перейдем в уравнениях (2) после преобразования
Коши к векторно-матричной форме:
+
(3)
Модель погрешностей ИНС можно описать следующим
векторным уравнением:
(4)
где X - вектор состояния, F - матрица динамики
системы , G - матрица влияния шумов системы, W - вектор белых шумов системы.
Получим вектор состояния X из системы (3):
Матрица динамики системы F, матрица влияния
шумов системы G и вектор белых шумов системы W примут следующий вид:
Для оценки влияния погрешностей проинтегрируем
выражение (4) методом Эллера:
Моделирование произведем в программном продукте
matlab. (Алгоритмы и листинги приведены в приложении)
После проведенных вычислений получены графики
изменения погрешностей координат и скоростей по времени (рис.1-6).
Рис. 3 - График погрешности долготного канала
Рис. 4 - График ошибки определения скорости в
долготном канале
Рис. 5 - График ошибки широтного канала
Рис. 6 - График ошибки определения скорости в
широтном канале
Рис. 7 - График ошибки определения высоты
Рис. 8 - График ошибки определения скорости в
высотном канале
Проанализировав графики можно увидеть, что на
графиках ошибок координаты и скорости долготного канала присутствует
расходимость. Она обусловлена влиянием нестабильного высотного канала через
перекрестные связи. Чтобы избежать данного негативного воздействия отделим
высотный канал, информация о котором будет поступать от внешних источников
(радиовысотомер).
В ходе принятых допущений преобразуем матрицу
динамики системы F обнулением элементов 6 строки, которая примет вид:
навигационный система алгоритм канал
В связи с принятыми изменениями внесем
соответствующие корректировки в алгоритм моделирования (приложение).
В итоге получим графики ошибок координат и
скоростей долготного и широтного каналов без воздействия на них вертикального
канала через перекрестные связи. Сравнив график ошибок долготного канала по
скорости и местоположению с воздействием через перекрестные связи высотного
канала и без воздействия высотного канала, можно увидеть что после изолирования
высотного канала пропала расходимость в долготном канале.
Рис. 9 - График погрешностей координат
долготного канала без воздействия высотного канала
Рис. 10 - График погрешностей определения
скорости в долготном канале без воздействия высотного канала
Рис. 11 - График определения координат широтного
канала без воздействия высотного канала
Рис. 12 - График определения скорости в широтном
канале без воздействия высотного канала
Выводы
После анализа уравнений ошибок и графиков
полученных в ходе моделирования можно увидеть, что при трех взаимосвязанных
каналах автономной ИНС нарастание погрешности со временем достигает величин,
при которых ИНС не будет отрабатывать достаточную точность. Данный процесс
фиксируется в расходимости графиков долготного канала (рис. 1-2). Это
обусловлено влиянием нестабильного вертикального канала через перекрестные связи,
которые можно проследить в матрице динамики системы F.
В итоге, чтобы исключить погрешности привносимые
вертикальным каналом в канал определения координат, высота не определяется в
вычислителе и исключены из общего алгоритма системы управления, используемые
для определения высоты и вертикальной скорости. Поскольку в алгоритме
вычисления координат и скорости необходимо использовать величины 
измерение
их осуществляется в радиовысотомере, и сигналы поступают в вычислитель ИНС.
После разрыва перекрестных связей вертикального канала и канала определения
координат результаты моделирования показывают отсутствие расходимости
долготного канала.
Используемая литература
1.
Навигационные приборы и системы (И.И. Помыкаев, В.П. Селезнев, Л.А.
Дмитроченко).
.
Лекции Антонова Д.А.
.
Ориентация и навигация подвижных объектов (Б.С. Алешин, К.К. Веремеенко, А.И.
Черноморский).
Приложение
[A,Y,G]=matr15() % функция инициализации матрицы
динамики системы A, вектора состояния Y, матрицы влияния шумов системы
% Входные параметры=0; % угол курса=0; % угол
широты=7.29e-5; % скорость вращения Земли=9.81; % ускорение свободного
падения=1.25e-3; % частота Шулера=0; % проекция абсолютной угловой скорости на
ось Х=u*cos(fi0); % проекция абсолютной угловой скорости на ось Y=u*sin(fi0); %
проекция абсолютной угловой скорости на ось Z=0; % производные =0; % проекций
абсолютных =0; % угловых скоростей на оси X,Y,Z=0; % ускорения измеренные=0; %
акселерометрами=-g0; % по осям X,Y,Z
% Начальные условия=15; % ошибка местоположения
по долготному каналу=0.1; % скоростная ошибка по долготному каналу=15; % ошибка
местоположения по широтному каналу=0.1; % скоростная ошибка по широтному
каналу=15; % ошибка местоположения по высотному каналу=0.1; % скоростная ошибка
по высотному каналу=3.48e-4; % Углы погрешности =3.48e-4; % построения
=3.48e-4; % базового 3-х гранника= 2.42e-8; % Проекции вектора
инструментальных= 2.42e-8; % и методических погрешностей измерителей угловой=
2.42e-8; % скорости на оси X,Y,Z= 0.005; % Погрешности измерения = 0.005; %
ускорения акселерометрами= 0.005; % по осям X,Y,Z
% Матрица динамики системы (для случая влияния через
перекрестные связи вертикального канала)
%A=[ 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0;
% (OMy^2+OMz^2-om0^2) 0 (dtOMz-OMx*OMy) 2*OMz
-(dtOMy+OMx*OMz) -2*OMy 0 -nz ny 1 0 0 0 0 0;
% 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0;
% -(dtOMz+OMx*OMy) -2*OMz (OMx^2+OMz^2-om0^2) 0
(dtOMx-OMy*OMz) 2*OMx nz 0 -nx 0 1 0 0 0 0;
% 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0;
% (dtOMy-OMx*OMz) 2*OMy -(dtOMx+OMy*OMz) 2*OMx
(2*om0^2+OMx^2+OMy^2) 0 -ny nx 0 0 0 1 0 0 0;
% 0 0 0 0 0 0 0 OMz -OMy 0 0 0 1 0 0;
% 0 0 0 0 0 0 -OMz 0 OMx 0 0 0 0 1 0;
% 0 0 0 0 0 0 OMy -OMx 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 0 0 0 0 0 0 0 0 0;
% 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
% 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
% 0 0 0 0 0 0 0 0 0 0 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 0;
(OMy^2+OMz^2-om0^2) 0 (dtOMz-OMx*OMy) 2*OMz
-(dtOMy+OMx*OMz) -2*OMy 0 -nz ny 1 0 0 0 0 0;
0 0 1 0 0 0 0 0 0 0 0 0 0 0;
(dtOMz+OMx*OMy) -2*OMz (OMx^2+OMz^2-om0^2) 0
(dtOMx-OMy*OMz) 2*OMx nz 0 -nx 0 1 0 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 1 0 0 0;
0 0 0 0 0 0 OMz -OMy 0 0 0 1 0 0;
0 0 0 0 0 -OMz 0 OMx 0 0 0 0 1 0;
0 0 0 0 0 OMy -OMx 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 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0;];
% Матрица динамики системы
%A = [0,1,0,0,0,0,0,0,0,0,0,0,0,0,0;...
%
1.3225*10^-10,0,0,0,0,-2.3*10^-5,0,10,0,1,0,0,0,0,0;...
% 0,0,0,1,0,0,0,0,0,0,0,0,0,0,0;...
% 0,0,-1.56*10^-6,0,0,0,-10,0,0,0,1,0,0,0,0;...
% 0,0,0,0,0,1,0,0,0,0,0,0,0,0,0;...
%
0,2.3*10^-5,0,0,3.125*10^-6,0,0,0,0,0,0,1,0,0,0;...
% 0,0,0,0,0,0,0,0,-1.15*10^-5,0,0,0,1,0,0;...
% 0,0,0,0,0,0,0,0,0,0,0,0,0,1,0;...
% 0,0,0,0,0,0,1.15*10^-5,0,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,0,0,0,0,0,0,0,0,0;...
% 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;...
% 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;...
% 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;...
% 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
% Вектор состояния
{x1,x2,x3,x4,x5,x6,a,b,g,dnx,dny,dnz,dOMx,dOMy,dOMz}
%Y =
%[15;0.1;15;0.1;15;0.1;0.005;0.005;0.005;0.009;0.009;0.009;0.009;0.009;0.009];
% Вектор
состояния=[x1;x2;x3;x4;x5;x6;a;b;g;dnx;dny;dnz;dOMx;dOMy;dOMz];
% Матрица влияния шумов системы=
[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,1,0,0,0;...
,0,0,1,0,0;...
,0,0,0,1,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,0,0,0,0];[tout, yout, Y]=eller(A, G, t0,
tfinal, y0, h) % функция реализации интегрирования методом Эллера =t0; y=y0;=t;
yout=y;(t<tfinal)=y+h*(A*y+G*wgn(6,1,20));=t+h;=[tout;t]; yout=[yout,y];=y;