Спутниковые системы навигации. Учебное пособие. Учебное пособие ( Лабораторный практикум на компьютере ) Київ 2008 1

Название	Учебное пособие ( Лабораторный практикум на компьютере ) Київ 2008 1
Дата	12.04.2023
Размер	4.49 Mb.
Формат файла
Имя файла	Спутниковые системы навигации. Учебное пособие.pdf
Тип	Учебное пособие #1057965
страница	3 из 14

1 2 3 4 5 6 7 8 9 ... 14

РАЗДЕЛ 2 Преобразование координат
2.1 Краткие сведения из теории
В спутниковых системах радионавигации применяются различные системы коорди- нат для расчета орбитального движения спутников и позиции потребителя. Системы ко- ординат, методы и алгоритмы их расчета и преобразования, на основании которых разра- ботаны программы данного раздела изложены в книге [1] (раздел 1. 3, стр. 33 -40), в руко- водстве [7]. Пакет программ в среде MatLab дается в папке COORDINATES и в прила- гаемых листингах программ.
Цель лабораторной работы: Изучение и практическое освоение систем координат, применяемых в спутниковых радионавигационных системах
2.2 Лабораторная работа 2. 1 «Преобразование координат»
Рекомендуется следующий порядок выполнения лабораторной работы.
1. Создайте папку COORDINATES_My и скопируйтев ее все программы из папки
COORDINATES.
2. Запустите MatLab.
3. Обратитесь к папке COORDINATES_My и откройте ее.
4. Откройте функцию ECEFLLH_N, внимательно изучите ее по комментариям и про- граммным процедурам, описывающим алгоритм расчета (раздел 1. 3 книги [1]).
5. Откройте пример расчета Pr_Coord1.m и выполните m- файл.
6. Основываясь на m- файле Pr_Coord1.m выполните задание 1.
7. Откройте функцию LLHECEF_N, внимательно изучите ее по комментариям и про- граммным процедурам, описывающим алгоритм расчета (раздел 1. 3 книги [1]).
8. Откройте пример расчета Pr_Coord2.m и выполните m- файл.
9. Основываясь на m- файле Pr_Coord2.m выполните задание 2.
10. Откройте функцию top_coord, внимательно изучите ее по комментариям и программ- ным процедурам, описывающим алгоритм расчета (раздел 1. 3 книги [1]).
11. Откройте пример расчета prim_top_coord.m и выполните m- файл.
12. Основываясь на m- файле prim_top_coord.m выполните задание 3.
13. Задание 1. По географической карте определите широту (градусы/минуты/секунды)и долготу (градусы/минуты/секунды) любого города Европы. Преобразуйте выбранные значения в радианы. Задайте высоту в метрах (произвольно, например, 195). В соот- ветствии с п. п. 5, 6 и выбранными входными данными сформируйте m-файл и выпол- ните его. Результат выполнения из командного окна MatLab перенесите в отчет.
38

14. Задание 2. Используя результаты, полученные в п. 13, в качестве входных данных в соответствии с п. п. 8, 9 сформируйте m-файл и выполните. Результат выполнения из командного окна MatLab перенесите в отчет.
15. Задание 3. Задайте координаты двух объектов, находящихся в прямой видимости ши- рота (градусы, минуты, секунды), долгота (градусы, минуты, секунды), высота (мет- ры). Преобразуйте заданные координаты широты и долготы в градусы. Используя эти исходные данные и п. п. 11, 12 сформируйте m-файл и выполните его. Результат вы- полнения из командного окна MatLab перенесите в отчет.
16. Обратитесь к папке .ECI_ECEF_LLH и откройте ее.
17. Откройте функцию eci_to_ecef, внимательно изучите ее по комментариям и программ- ным процедурам, описывающим алгоритм расчета (раздел 1. 3 книги [1]).
18. Откройте пример расчета Pr_eci_ecef.m и выполните m- файл.
19. Основываясь на m- файле Pr_eci_ecef.m выполните задание 4.
20. Откройте функцию llh_to_eci, внимательно изучите ее по комментариям и программ- ным процедурам, описывающим алгоритм расчета (раздел 1. 3 книги [1]).
21. Откройте пример расчета Pr_ecef_eci.m и выполните m- файл.
22. Основываясь на m- файле Pr_ecef_eci.m выполните задание 5.
23. Задание 4. Задайтекоординаты и скорости объекта satpos_eci, истинное звездное вре- мя s0, текущее время ti. Используя эти исходные данные и п. п. 17, 18 сформируйте m- файл и выполните его. Результат выполнения из командного окна MatLab перенесите в отчет.
24. Задание 5. Выберите эллипсоид и задайте его полуоси a, b; задайте текущее время
ti,истинное звездное время time_s0, координаты потребителя llh_loc. Используя эти исходные данные и п. п. 20, 21 сформируйте m-файл и выполните его. Результат вы- полнения из командного окна MatLab перенесите в отчет.
25. Обратитесь к папке . TEST и откройте ее.
26. Откройте функции ECEF_to_LLH_Dg_Zu, ECEF_to_LLH_Itera, ECEF_to_LLH_Kelly, внимательно изучите их по комментариям и программным процедурам, описывающим алгоритм расчета (раздел 1. 3 книги [ ], руководство [ ]).
27. Откройте пример тестирования расчетов Test_Coord.m , изучите его, выполните m- файл.
28. Основываясь на m- файле Test_Coord.m выполните задание 6.
29. Задание 6. Используя входные данные из заданий 1- 5 и функции
ECEF_to_LLH_Dg_Zu, ECEF_to_LLH_Itera, ECEF_to_LLH_Kelly выполните сопостав- ление расчетов по точной , итерационной и приближенной формулам при изменение
39

высоты в пределах 0 – 500 м; 1 – 10 км; 1 – 20 000 км. Расчеты проведите в области эк- ватора, северного полюса и Киева. Результаты выполнения из графиков включите в отчет.
2.3 Вопросы и задания для самоподготовки
1. Какие системы координат применяются в спутниковых радионавигационных систе- мах?
2. Какая разница между геоцентрическими и геодезическими координатами?
3. Что обозначают понятия правая и левая системы координат?
4. Что обозначают понятия подвижная и неподвижная системы координат?
5. Запишите в аналитическом виде формулы перехода из пространственной эллипсоид- ной географической системы в геоцентрическую фиксированную систему (ECEF).
6. Даете определение пространственной эллипсоидной географической системе коорди- нат (центр, широта, долгота, высота).
7. Как определяются эллипсоид, геоид?
8. Сформулируйте определение системы координат, имеющей международное обозначе- ние ECEF.
9. Сформулируйте определение системе координат, имеющей международное обозначе- ние ECI.
10. Дайте определение топоцентрической системе координат (цент, направления осей).
11. Дайте определение системе координат WGS 84 (цент, направления осей, параметры эллипсоида, в каких спутниковых радионавигационных системах является опорной).
12. Дайте определение системе координат ПЗ 90 (цент, направления осей, параметры эл- липсоида, в каких спутниковых радионавигационных системах является опорной).
13. Объясните понятие «прямая видимость».
14. Запишите формулу перевода градусов, минут, секунд в градусы, радианы; составьте программу в виде m- файла и убедитесь в правильности работы программы.
15. Запишите формулу перевода радиан в градусы, минуты, секунды; в градусы. Составьте программу в виде m- файла и убедитесь в правильности работы программы.
40

2.4 Тексты программ
2.4.1 Функции и файлы из папки COORDINATES
Функция ECEFLLH_N
function [XYZ] = ECEFLLH_N(llh,ab)
%Имя функции: ECEFLLH_N
%Назначение функции: преобразование координат из географической системы в прямоугольную
%Входные данные:
%llh.lon-долгота;
%llh.lat-широта;
%llh.h-высота;
%ab.a-большая полуось эллипсоида;
%ab.b- малая полуось эллипсоида в WGS-84;
%Выходные данные:
%XYZ.x,XYZ.y,XYZ.z- координаты X, Y, Z соответственно в ECEF
% Справочные данные:
%ECEF- прямоугольная геоцентрическая система координат
%a=6378137.0 (м)- большая полуось эллипсоида для WGS-84;
%b=6356752.314 (м)- малая полуось эллипсоида для WGS-84;
%A_PZ90_M =6 378 136 (м)- большая полуось эллипсоида для ПЗ 90;
%B_PZ90_M = 6356751.36174 (м)- малая полуось эллипсоида для ПЗ 90; a2=ab.a*ab.a; b2=ab.b*ab.b; r=a2/sqrt(a2*cos(llh.lat)*cos(llh.lat)+b2*sin(llh.lat)*sin(llh.lat));
XYZ.x=(r+llh.h)*cos(llh.lat)*cos(llh.lon);
XYZ.y=(r+llh.h)*cos(llh.lat)*sin(llh.lon);
XYZ.z=(b2/a2*r+llh.h)*sin(llh.lat);
Файл Pr_Coord1.m
%Имя m- файла:Pr_Coord1.m
%Пример расчета llh.lat=0.881278698506528;llh.lon=0.53169758803674; llh.h=122.899802776054; ab.a=6378137.0; ab.b=6356752.314;
[XYZ] = ECEFLLH_N(llh,ab)
Функция LLHECEF_N
function [llh] = LLHECEF_N(XYZ,ab)
%Имя функции: LLHECEF_N
%Назначение функции: преобразование координат из прямоугольной системы в географическую
%Входные данные:
41

%XYZ.x,XYZ.y,XYZ.z- координаты X, Y, Z соответственно в ECEF
%ab.a-большая полуось эллипсоида;
%ab.b- малая полуось эллипсоида в WGS-84;
%Выходные данные:
%llh.lon-долгота;
%llh.lat-широта;
%llh.h-высота;
% Справочные данные:
%ECEF- прямоугольная геоцентрическая система координат
%a=6378137.0 (м)- большая полуось эллипсоида для WGS-84;
%b=6356752.314 (м)- малая полуось эллипсоида для WGS-84;
%A_PZ90_M =6 378 136 (м)- большая полуось эллипсоида для ПЗ 90;
%B_PZ90_M = 6356751.36174 (м)- малая полуось эллипсоида для ПЗ 90; a=6378137.0; b=6356752.314; a2=ab.a*ab.a; b2=ab.b*ab.b; xy = sqrt(XYZ.x*XYZ.x + XYZ.y*XYZ.y); thet = atan(XYZ.z*ab.a/(xy*ab.b)); esq = 1.0-b2/a2; epsq = a2/b2-1.0; llh.lat = atan((XYZ.z+epsq*ab.b*(sin(thet)^3))/(xy-esq*ab.a*(cos(thet)^3))); llh.lon = atan2(XYZ.y,XYZ.x);%! if llh.lon < 0 llh.lon = 2*pi + llh.lon; end ; r = a2/sqrt(a2*cos(llh.lat)*cos(llh.lat) + b2*sin(llh.lat)*sin(llh.lat)); llh.h = xy/cos(llh.lat)-r; end
Файл :Pr_Coord2.m
%Имя m- файла:Pr_Coord2.m
%Пример расчета ab.a=6378137.0; ab.b=6356752.314;
XYZ.x=3.504451023000798e+006;XYZ.y=2.061316876000462e+006;
XYZ.z=4.897990974997338e+006;
[llh] = LLHECEF_N(XYZ,ab)
Функция top_coord
function [top] = top_coord(rec_llh, rec_xyz, nlo_xyz)
42

% Имя функции: top_coord
%Назначение функции: расчет топоцентрических координат объекта по заданным
%географическим (долгота, широта, высота) и геоцентрическим (x, y, z)
%координатам приемника и геоцентрическим координатам объекта (x, y, z)
% Входные данные:
% rec_llh.lat - широта (рад) приемника;
%rec_llh.lon -- долгота (рад) приемника;
%rec_llh.h- высота (м) приемника;
%прямоугольные геоцентрические координаты приемника (м):
% rec_xyz.x
% rec_xyz.y
%,, rec_xyz.z
%прямоугольные геоцентрические координаты объекта (м):
% ns.x - координата x;
% ns.y -координата y;
% ns.z- координата z ;
% Выходные данные:
% top.s - проекция вектора дальности на ось (м) , направленную на Юг (South)
% top.e - проекция вектора дальности на ось (м) , направленную на Восток (East)
% top.z - проекция вектора дальности на ось (м) , направленную в Зенит
% top.daln - дальность до объекта (м)
% top.az - угол азимута объекта (градус)
% top.el - угол видимости объекта (градус) rx = nlo_xyz.x - rec_xyz.x; ry = nlo_xyz.y - rec_xyz.y; rz = nlo_xyz.z - rec_xyz.z; r_sat = sqrt(rx*rx + ry*ry + rz*rz); r_rec = sqrt((rec_xyz.x)^2 + (rec_xyz.y)^2+ (rec_xyz.z)^2); top.r = r_sat; rx1 = rx; ry1 = ry; rz1 = rz; sin_lat = sin(rec_llh.lat); cos_lat = cos(rec_llh.lat); sin_lon = sin(rec_llh.lon); cos_lon = cos(rec_llh.lon);
% Projections of vector of range in topocentric coordinate system: top.e = -sin_lon * rx1 + cos_lon * ry1; top.s = cos_lon * sin_lat * rx1 + sin_lon * sin_lat * ry1 - cos_lat * rz1; top.z = cos_lat * cos_lon * rx1 + cos_lat * sin_lon * ry1 + sin_lat * rz1;
% azimut: отсчет по часовой стрелке от оси направленной на Север (N or -S) (-top.s) eps = 10e-10; if ( (abs(top.e) < eps) || (abs(top.s) < eps))
43

top.az = 0.0; else top.az = atan2(top.e,-top.s); end; if (top.az < 0.0) top.az = top.az + pi * 2; end;
% elevation: cos_el_top = (rec_xyz.x * rx + rec_xyz.y * ry + rec_xyz.z * rz) / (r_sat * r_rec); if ( cos_el_top >= 1.00 ) el = 0.0; else if ( cos_el_top <= -1.00 ) el = pi; else el = acos(cos_el_top); end; end; top.el = pi / 2.0 - el;
Файл prim_top_coord.m
%Имя m-файла: prim_top_coord.m
%Пример расчета a=6378137.0; b=6356752.314; % для WGS-84;
% Коэффициенты перевода градусов в радианы и обратно
A2R = pi/180;
R2A = 180/pi;
%Входные данные координаты, например, приемника rec_deg.lon = 100; rec_deg.lat = 40; rec_deg.h = 0;
%Входные данные координат объекта nlo_deg.lon = 280; nlo_deg.lat = -40; nlo_deg.h = 0;
%Преобразование градусов в радианы rec_llh.lon = rec_deg.lon * A2R; rec_llh.lat = rec_deg.lat * A2R; rec_llh.h = rec_deg.h; nlo_llh.lon = nlo_deg.lon * A2R; nlo_llh.lat = nlo_deg.lat * A2R; nlo_llh.h = nlo_deg.h;
44

%Преобразование координат приемника и объекта систему ECEF
[rec_xyz] = ECEFLLH(a, b, rec_llh);
[nlo_xyz] = ECEFLLH(a, b, nlo_llh);
%Преобразование координат приемника и объекта в топоцентрическую
%систему координат
[top] = top_coord(rec_llh, rec_xyz, nlo_xyz);
%Вывод данных приемника в топоцентрической системе координат fprintf('e=%22.16e s=%22.16f z=%22.16f az=%f el=%f r=%f \n', top.e, top.s, top.z, top.az*R2A, top.el*R2A, top.r);
%Вывод данных объекта в топоцентрической системе координат fprintf('e=%22.16e s=%22.16f z=%22.16f az=%f el=%f r=%f \n', top.e, top.s, top.z, top.az*R2A, top.el*R2A, top.r);
Функция ECEFLLH
function [R] = ECEFLLH(a, b, llh)
%Имя функции: ECEFLLH
%Назначение- вариант функции ECEFLLH_N a2 = a * a; b2 = b * b; n = a2 / sqrt(a2 * cos(llh.lat)*cos(llh.lat) + b2 * sin(llh.lat) * sin(llh.lat));
R.x = (n + llh.h) * cos(llh.lat) * cos(llh.lon);
R.y = (n + llh.h) * cos(llh.lat) * sin(llh.lon);
R.z = (b2 / a2 * n + llh.h) * sin(llh.lat);
2.4.2 Функции и файлы из папки ECI_ECEF_LLH
Функция eci_to_ecef
function [satpos_ecef] =eci_to_ecef(s0, ti, satpos_eci)
%Имя функции:eci_to_ecef
%Функция преобразования координат
%Входные данные: s0 - истинное звездное время в текущий момент обсервации ,
%ti - текущее время; satpos_eci
%Структура satpos_eci
%satpos_eci.x - координата x в абсолютной неподвижной системе координат (ECI);
%satpos_eci.y - координата y в абсолютной неподвижной системе координат (ECI);
%satpos_eci.z - координата z в абсолютной неподвижной системе координат (ECI);
%satpos_eci.vx - скорость vx в абсолютной неподвижной системе координат (ECI);
%satpos_eci.vy - скорость vy в абсолютной неподвижной системе координат (ECI);
% satpos_eci.vz - скорость vz в абсолютной неподвижной системе координат (ECI);
%Выходные данные:
% Структура satpos_ecef
%satpos_ecef.x - координата x в подвижной системе координат (ECEF);
45

%satpos_ecef.y - координата y в подвижной системе координат (ECEF);
%satpos_ecef.z - координата z в подвижной системе координат (ECEF);
%satpos_ecef.vx - скорость по оси x в подвижной системе координат (ECEF);
%satpos_ecef.vy - скорость по оси z в подвижной системе координат (ECEF);
%satpos_ecef.vz- скорость по оси z в подвижной системе координат (ECEF);
%Коэффициенты
% SEC_IN_RAD - коэффициент преобразования секунд в радианы
% s0(radian) = s0 (sek) * SEC_IN_RAD, where
% SEC_IN_RAD = 2 * pi / (24 * 3600) = pi / 43200
SEC_IN_RAD = pi / 43200;
OMEGA_Z = 0.7292115e-4; %( скорость вращения Земли (angular speed of rotation of the Earth, рад/cek) s_zv = s0 * SEC_IN_RAD + OMEGA_Z * ti; cos_s = cos(s_zv); sin_s = sin(s_zv); satpos_ecef.x = satpos_eci.x * cos_s + satpos_eci.y * sin_s; satpos_ecef.y = -satpos_eci.x * sin_s + satpos_eci.y * cos_s; satpos_ecef.z = satpos_eci.z; satpos_ecef.vx = satpos_eci.vx * cos_s + satpos_eci.vy * sin_s + OMEGA_Z * satpos_ecef.y; satpos_ecef.vy = -satpos_eci.vx * sin_s + satpos_eci.vy * cos_s - OMEGA_Z * satpos_ecef.x; satpos_ecef.vz = satpos_eci.vz;
Файл Pr_eci_ecef.m
%Имя файла:Pr_eci_ecef.m
%Назначение- пример преобразования координат
%Входные данные: s0 = 400; ti =500; satpos_eci.x= 20000000;% координата x в абсолютной неподвижной системе координат (ECI); satpos_eci.y= 15000000 ;%координата y в абсолютной неподвижной системе координат (ECI); satpos_eci.z = 10000000;% координата z в абсолютной неподвижной системе координат (ECI); satpos_eci.vx = 5000;% скорость vx в абсолютной неподвижной системе координат (ECI); satpos_eci.vy= 6000;% скорость vy в абсолютной неподвижной системе координат (ECI); satpos_eci.vz= 7000;%скорость vz в абсолютной неподвижной системе координат (ECI);
%Выходные данные
[satpos_ecef] =eci_to_ecef(s0, ti, satpos_eci)
Функция llh_to_eci
function [eci_llh, eci_xyz] = llh_to_eci(a, b, ti, time_s0, llh_loc) ;
%Имя функции:llh_to_eci
46

%Функция вычисляет позицию приемника в абсолютной геоцентрической системе координат (ECI)
%Входные данные:
%a, b-большая и малая полуоси земного эллипсоида (метр);
%ti- текущее время (секунды) ,
%time_s0- истинное звездное время,
%Структура llh_loc - координаты приемника; {
%llh_loc.lon-долгота (радиан);
%llh_loc.lat-широта (радиан);
%llh_loc.h-высота (метр);
%Выходные данные:
%Структура eci_llh - географические координаты приемника в абсолютной геоцентрической системе координат (ECI)
%eci_llh.lon - долгота (радиан);
%eci_llh.lat - широта (радиан);
%eci_llh.h - высота (метр);
%Структура eci_xyz- координаты приемника в абсолютной прямоугольной геоцентрической системе координат (ECI)
% eci_xyz.x - координата x;
% eci_xyz.y - координата y;
% eci_xyz.z - координата z;
OMEGA_Z = 0.7292115e-4; % угловая скорость вращения Земли, (рад/cek)
SEC_IN_RAD = 7.2722052166430e-005; % (PI / 43200.0 ) // Number radian eci_llh.lon = llh_loc.lon + ti * OMEGA_Z + SEC_IN_RAD * time_s0; eci_llh.lat = llh_loc.lat; eci_llh.h = llh_loc.h; eci_xyz = llh_to_ecef( a, b, eci_llh);
Файл :Pr_ecef_eci.m
%Имя файла:Pr_ecef_eci.m
%Назначение- пример преобразования координат
%Входные данные: a=6378137.0;% (м)- большая полуось эллипсоида для WGS-84; b=6356752.314;% (м)- малая полуось эллипсоида для WGS-84; ti= 700;% текущее время (секунды) , time_s0= 100;% истинное звездное время,
%координаты приемника: llh_loc.lon = 30*pi/180; %долгота (радиан); llh_loc.lat= 55*pi/180;%широта (радиан); llh_loc.h= 184; %высота (метр);
%Выходные данные
47

[eci_llh, eci_xyz] = llh_to_eci(a, b, ti, time_s0, llh_loc)
Функция llh_to_ecef
function [XYZ] = llh_to_ecef( a, b, llh)
%Имя функции:LLH_to_ECEF.m
%Функция преобразования географических координат в прямоугольную геоцентрическую систему координат (ECEF)
%Входные данные:
%Структура llh
%llh.lon-долгота (радиан),
%llh.lat-широта (радиан),
%llh.h-высота (метр);
%a, b-большая и малая полуоси земного эллипсоида в WGS-84 (метр);
%Выходные данные:
%Структура XYZ
%XYZ.x - координата x в ECEF;
%XYZ.y - координата y в ECEF;
%XYZ.z - координата z в ECEF; a2=a * a; b2=b * b; r = a2 / sqrt(a2 * cos(llh.lat) * cos(llh.lat) + b2 * sin(llh.lat) * sin(llh.lat));
XYZ.x = (r + llh.h) * cos(llh.lat) * cos(llh.lon);
XYZ.y = (r + llh.h) * cos(llh.lat) * sin(llh.lon);
XYZ.z = (b2 / a2 * r + llh.h) * sin(llh.lat);
2.4.3 Функции и файлы из папки TEST
Функция ECEF_to_LLH_Dg_Zu
function [llh_D] = ECEF_to_LLH_Dg_Zu(a, b, XYZ)
%Имя функции: ECEF_to_LLH_Dg_Zu
%Назначение функции преобразование координат из прямоугольной системы в географическую
% по точным формулам
%Входные данные:
%XYZ.x- координата x, (м) в системе ECEF;
%XYZ.y- координата y, (м) в системе ECEF;
%XYZ.z- координата z, (м) в системе ECEF;
%a=6378137.0 (м)- большая полуось эллипсоида для WGS-84;
%b=6356752.314 (м)- малая полуось эллипсоида для WGS-84;
%A_PZ90_M =6 378 136 (м)- большая полуось эллипсоида для ПЗ 90;
%B_PZ90_M = 6356751.36174 (м)- малая полуось эллипсоида для ПЗ 90;
%Выходные данные:
48

%llh_D.lambda-долгота;
%llh_D.Fi-широта;
%llh_D.h-высота; e2=0.00669437999; %квадрат эксцентриситета эллипсоида a=6378137;% экваториальный радиус эллипсоида b=a*sqrt(1-e2);% полярный радиус эллипсоида
%{
XYZ.x= 1.449310528799138e+007;
XYZ.y= 8.513116350113131e+006;
XYZ.z= 2.031312580583197e+007;
%} w=sqrt(XYZ.x*XYZ.x+XYZ.y*XYZ.y); l=e2/2; m=(w/a)*(w/a); n= ((1-e2)*XYZ.z/b)^2; i=-(2*l*l+m+n)/2; k= l*l*(l*l-m-n); q=(m+n-4*l*l)^3/216 +m*n*l*l;
D=sqrt((2*q-m*n*l*l)*m*n*l*l); bet=i/3-(q+D)^(1/3) - (q-D)^(1/3); t=sqrt(sqrt(bet*bet-k)-(bet+i)/2)-sign(m-n)*sqrt((bet-i)/2); w1=w/(t+l); z1=(1-e2)*XYZ.z/(t-l); llh_D.Fi= atan(z1/((1-e2)*w1)); llh_D.lambda= 2*atan((w-XYZ.x)/XYZ.y); llh_D.h= sign(t-1+l)*sqrt((w-w1)^2 +(XYZ.z-z1)^2);
Функция ECEF_to_LLH_Itera
function [llh] = ECEF_to_LLH_Itera(a, b, XYZ)
%Имя функции: ECEF_to_LLH_Itera
%Назначение функции преобразование координат из прямоугольной системы в географическую
% по итерационным формулам
%Входные данные:
%XYZ.x- координата x, (м) в системе ECEF;
%XYZ.y- координата y, (м) в системе ECEF;
%XYZ.z- координата z, (м) в системе ECEF;
%a=6378137.0 (м)- большая полуось эллипсоида для WGS-84;
%b=6356752.314 (м)- малая полуось эллипсоида для WGS-84;
%A_PZ90_M =6 378 136 (м)- большая полуось эллипсоида для ПЗ 90;
%B_PZ90_M = 6356751.36174 (м)- малая полуось эллипсоида для ПЗ 90;
%Выходные данные:
%llh.lon - долгота;
49

% llh.lat - широта;
%llh.h- высота; a2 = a * a; b2 = b * b; sqrt_xy = sqrt(XYZ.x * XYZ.x + XYZ.y*XYZ.y); e2 = 1.0 - b2 / a2; % ' e =Эксцентриситет, e2 = e^2 eps = 0.001; delta_h = 100; n = 0; v = a;% радиус кривизны в главном вертикале h = 0; while (abs(delta_h) > eps) n = n + 1; fi = atan(XYZ.z / (sqrt_xy * (1 - e2 * v / (v + h)))); % lat llh.lat = fi; sin_fi = sin(fi); cos_fi = cos(fi); v = a / sqrt(1 - e2 * sin_fi * sin_fi); llh.lon = atan2(XYZ.y,XYZ.x); llh.h = sqrt_xy / cos_fi - v; delta_h = llh.h - h;
% fprintf('n=%i h=%f lat=%f lon=%f h=%f delta_h= %f \n',n, h, llh.lat, llh.lon, llh.h, delta_h); h = llh.h; end; end
Функция ECEF_to_LLH_Kelly
function [llh] = ECEF_to_LLH_Kelly(a, b, XYZ)
%Имя функции: ECEF_to_LLH_Kelly
%Назначение функции преобразование координат из прямоугольной системы в географическую
% по итерационным формулам
%Входные данные:
%XYZ.x- координата x, (м) в системе ECEF;
%XYZ.y- координата y, (м) в системе ECEF;
%XYZ.z- координата z, (м) в системе ECEF;
%a=6378137.0 (м)- большая полуось эллипсоида для WGS-84;
%b=6356752.314 (м)- малая полуось эллипсоида для WGS-84;
%A_PZ90_M =6 378 136 (м)- большая полуось эллипсоида для ПЗ 90;
%B_PZ90_M = 6356751.36174 (м)- малая полуось эллипсоида для ПЗ 90;
%Выходные данные:
%llh.lon - долгота;
% llh.lat - широта;
50

%llh.h- высота; a2 = a * a; b2 = b * b; xy = sqrt(XYZ.x * XYZ.x + XYZ.y*XYZ.y); thet = atan(XYZ.z * a / (xy * b)); esq = 1.0 - b2/a2; % ' e =Эксцентриситет, esq = e^2 epsq = a2/b2 - 1.0; llh.lat = atan((XYZ.z + epsq * b * (sin(thet)^3)) / (xy - esq * a * (cos(thet)^3))); llh.lon = atan2(XYZ.y,XYZ.x); if llh.lon < 0 llh.lon = 2*pi + llh.lon; end ; r = a2 / sqrt(a2 * cos(llh.lat) * cos(llh.lat) + b2 * sin(llh.lat) * sin(llh.lat)); llh.h = xy/cos(llh.lat) - r; end
Функция ECEFLLH
function [R] = ECEFLLH(a, b, llh)
%Имя функции: ECEFLLH
%Назначение- вариант функции ECEFLLH_N a2 = a * a; b2 = b * b; n = a2 / sqrt(a2 * cos(llh.lat)*cos(llh.lat) + b2 * sin(llh.lat) * sin(llh.lat));
R.x = (n + llh.h) * cos(llh.lat) * cos(llh.lon);
R.y = (n + llh.h) * cos(llh.lat) * sin(llh.lon);
R.z = (b2 / a2 * n + llh.h) * sin(llh.lat);
Файл Test_Coord.m
%Имя m-файла:Test_Coord.m
%Назначение: пример тестирования программ преобразования координат
%llh.lat = 0.88032730015257; %50 град; 26 мин.; 20.54 с
%ввод входных данных llh.lat =55*pi/180; llh.lon = 0.53109641675259;%30 град; 25 мин.; 46.4995 с a=6378137.0; b=6356752.314;
% ввод высоты и шага изменения высоты llh0.h=0;%184;%высота в метрах step_h= 1000; kol=40000;
%nn=0;
51

for nn=1:kol
% nn=nn+1; llh.h= llh0.h+step_h*(nn-1);
%применение функции вычисления прямоугольных координат при переменной высоте
[R] = ECEFLLH(a, b, llh);
%преобразование прямоугольных координат в географические по точной формуле
[llh_D] = ECEF_to_LLH_Dg_Zu(a, b, R);
%преобразование прямоугольных координат в географические поприближенной формуле
%указание снять блокировку строки и переписать вывод в графике
% [llh] = ECEF_to_LLH_Kelly(a, b,R);
%[llh_K] =Kelly(a, b, R);
%преобразование прямоугольных координат в географические по итерационной формуле
[llh1] = ECEF_to_LLH_Itera(a, b, R);
%{ delta.lat(nn)=llh_D.Fi- llh.lat; delta.lon(nn)=llh_D.lambda- llh.lon; delta.h(nn)= llh_D.h- llh.h;
%} h1(nn)= llh.h; %llh0.h+step_h*(nn-1); delta.lat(nn) = llh1.lat - llh.lat;
%delta.lon(nn)=llh_D.lambda- llh_K.lon;
%delta.h(nn) = llh.h - llh_K.h; delta.h(nn) = llh.h - llh1.h; end
%Графика subplot (1,2,1), plot(h1(2:kol)/1000,delta.h(2:kol)*1000), grid on xlabel('a','FontSize',14,'FontName','TimesNewRoman') set(get(gcf,'CurrentAxes'),'FontSize',14,'FontName','TimesNewRoman'); subplot (1,2,2),,plot(h1(2:kol)/1000,delta.lat(2:kol)*180*3600/pi), grid on xlabel('б','FontSize',14,'FontName','TimesNewRoman') set(get(gcf,'CurrentAxes'),'FontSize',14,'FontName','TimesNewRoman');
%plot(h1,delta.lat*180/pi)
%plot(h1*10^(-3),delta.h)
%plot(1:kol+1, h1)
52

РАЗДЕЛ 3 Время
3.1 Краткие сведения из теории
В спутниковой радионавигации время имеет большое значение, поскольку основные навигационные определения производятся по формулам, в которых параметр времени присутствует многократно. Прежде всего, это время распространения электромагнитного сигнала от навигационного спутника до потребителя, время «включения» часов спутника, время синхронизации данных передаваемых со спутника, время прохождения электромаг- нитного сигнала через атмосферу, влияние на время релятивистских эффектов, совмеще- ние шкал времени спутника и потребителя, системное время СРНС, опорные моменты от- счета времени (эпохи), единицы счета времени (год, неделя, день, час, минута, секунда, миллисекунда) и многое другое.
Алгоритмы расчета времени, запрограммированные в прилагаемом пакете программ, изложены в книге [1] (раздел 1. 3, стр.40 -50), а также в источниках, приведенных в биб- лиографии к книге [1]. Пакет программ в среде MatLab дается в папке TIME и в прила- гаемых листингах программ.
Цель лабораторной работы: Изучение основных временных составляющих, при- меняемых в алгоритмах и программах спутниковой аппаратуры потребителя для решений навигационных задач.
3.2 Лабораторная работа 3. 1 «Время в спутниковых радионавигационных сис-

1 2 3 4 5 6 7 8 9 ... 14