Разработка и исследование характеристик платформенной инерциальной навигационной системы полуаналитического типа

Цель настоящей работы разработать алгоритм платформенной инерциальной навигационной системы, работающей в геоцентрической системе координат, и определяющей в этой системе следующие параметры:

Координаты Скорости Углы ориентации Так же предусмотрена задача получения позиционной и скоростной информации в географической системе координат (ц, л,), вычисление углов курса, крена и тангажа.

Исследование точностных характеристик системы по уравнениям ошибок и оценка влияния ошибок начальной выставки и гироскопов на точность ИНС производится при помощи алгоритма разработанного в программной среде matlab.

В конструкции ИНС используется акселерометры А-17 и лазерные гироскопы ГЛ-1 производства Раменского приборостроительного завода.

Исходные данные

Траекторные условия: полет с постоянной скоростью W=900км/ч, на постоянной высоте H=10 000м, курс постоянный произвольный, время полета 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 означает, что

R = 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, т. е. на величину разности географической и геоцентрической широты. Взаимное расположение трехгранников определяется таблицей направляющих косинусов:

x2 y2 z2

x11 0 0;

y10 cos (ф ф1) -sin (ф ф1);(1.20)

z10 sin (ф ф1) cos (ф ф1).

Выражение для разности (ф ф1) определится:

Вследствие малости величин и, считая также величину h/a малой и раскладывая правые части указанных формул в ряды по степеням и h/a, будем иметь Ввиду малости и упрощается матрица направляющих косинусов. Принимая cos () = 1, sin) = получим

x2 y2 z2

x11 0 0;

y10 1 ;

z10 1

При подстановке значения = 0,0067 получаем максимальное отклонение истинной вертикали от геоцентрической, равное = 0,335, что соответствует и имеет место на широте. С увеличением 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. Навигационные приборы и системы (И.И. Помыкаев, В. П. Селезнев, Л.А. Дмитроченко).

2. Лекции Антонова Д.А.

3. Ориентация и навигация подвижных объектов (Б.С. Алешин, К. К. Веремеенко, А.И. Черноморский).

Приложение

function [A, Y, G]=matr15() % функция инициализации матрицы динамики системы A, вектора состояния Y, матрицы влияния шумов системы

% Входные параметры

psi0=0; % угол курса

fi0=0; % угол широты

u=7.29e-5; % скорость вращения Земли

g0=9.81; % ускорение свободного падения

om0=1.25e-3; % частота Шулера

OMx=0; % проекция абсолютной угловой скорости на ось Х

OMy=u*cos (fi0); % проекция абсолютной угловой скорости на ось Y

OMz=u*sin (fi0); % проекция абсолютной угловой скорости на ось Z

dtOMx=0; % производные

dtOMy=0; % проекций абсолютных

dtOMz=0; % угловых скоростей на оси X, Y, Z

nx=0; % ускорения измеренные

ny=0; % акселерометрами

nz=-g0; % по осям X, Y, Z

% Начальные условия

x1=15; % ошибка местоположения по долготному каналу

x2=0.1; % скоростная ошибка по долготному каналу

x3=15; % ошибка местоположения по широтному каналу

x4=0.1; % скоростная ошибка по широтному каналу

x5=15; % ошибка местоположения по высотному каналу

x6=0.1; % скоростная ошибка по высотному каналу

a=3.48e-4; % Углы погрешности

b=3.48e-4; % построения

g=3.48e-4; % базового 3-х гранника

dOMx= 2.42e-8; % Проекции вектора инструментальных

dOMy= 2.42e-8; % и методических погрешностей измерителей угловой

dOMz= 2.42e-8; % скорости на оси X, Y, Z

dnx = 0.005; % Погрешности измерения

dny = 0.005; % ускорения акселерометрами

dnz = 0.005; % по осям X, Y, Z

% Матрица динамики системы (для случая влияния через перекрестные связи вертикального канала)

%A=[ 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0;

% (OMy²+OMz²-om0²) 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²+OMz²-om0²) 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²+OMx²+OMy²) 0 -ny nx 0 0 0 1 0 0 0;

% 0 0 0 0 0 0 0 OMzOMy 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 OMyOMx 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;];

% Матрица динамики системы (для случая отсутствия влияния вертикального канала)

A=[ 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0;

(OMy²+OMz²-om0²) 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²+OMz²-om0²) 0 (dtOMx-OMy*OMz) 2*OMx nz 0 -nx 0 1 0 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 1 0 0 0;

0 0 0 0 0 0 0 OMzOMy 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 OMyOMx 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;];

% Матрица динамики системы

%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];

% Вектор состояния

Y=[x1;x2;x3;x4;x5;x6;a;b;g;dnx;dny;dnz;dOMx;dOMy;dOMz];

% Матрица влияния шумов системы

G = [0,0,0,0,0,0;…

1,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,1,0,0,0;…

0,0,0,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,0,0,0;…

0,0,0,0,0,0;…

0,0,0,0,0,0];

function [tout, yout, Y]=eller (A, G, t0, tfinal, y0, h) % функция реализации интегрирования методом Эллера

t=t0; y=y0;

tout=t; yout=y;

while (t

y=y+h*(A*y+G*wgn (6,1,20));

t=t+h;

tout=[tout;t]; yout=[yout, y];

end

Y=y;