Спутниковые системы навигации. Учебное пособие. Учебное пособие ( Лабораторный практикум на компьютере ) Київ 2008 1
 Скачать 4.49 Mb.
|
|
4.2 Размножения эфемерид спутников ГЛОНАСС (иллюстрация решения сис- темы дифференциальных уравнений) 4.2.1 Краткие сведения из теории Одним из отличий системы ГЛОНАСС от GPS является то, что размножение эфеме- рид в системе ГЛОНАСС производится через решение системы дифференциальных урав- нений. В полном виде система дифференциальных уравнений для ГЛОНАСС приводится в интерфейсном контрольном документе [ ], упрощенная принципиальная схема решения дается в книге [1] ( параграфы 1. 3. 4, стр. 51- 56, 4. 2.1, стр. 198-201 ). Цель лабораторной работы: Овладение методом решения дифференциальных уравнения для размножения эфемерид спутниковой системы ГЛОНАСС. 4.2.2 Лабораторная работа 4. 3 «Решения системы дифференциальных уравнений» Приведенная программа иллюстрирует применение функции MatLab [7, 8] при ре- шении системы дифференциальных уравнений методом Рунге- Кутта для размножения эфемерид спутников ГЛОНАСС. На CD- диске программа расположена в папке ОРБИ- ТА_GLONASS. Для выполнения работы в качестве входных данных потребуются коор- динаты и скорости спутников ГЛОНАСС, которые являются начальными условиями при 96 решении системы дифференциальных уравнений; время, на которое заданаются эти пара- метры в данных, передаваемых со спутников, время на которое рассчитываются коорди- наты и скорости навигационных спутников. Поскольку приводится существенно упро- щенный алгоритм размножения эфемерид, то другие составляющие данных с навигацион- ного спутника не учитываются. Рекомендуется следующий порядок выполнения лабораторной работы. 1. Создайте папку ОРБИТА_GLONASS_My и скопируйте в ее из папки ОРБИ- ТА_GLONASS m-файлОRBITA_1.mи функцию orbit_GL. 2. Изучите программные процедуры и комментарии к m-файлуОRBITA_1.mи функции orbit_GL. Выполните задание 1. 3. Задание 1. С сайта Российского космического агентства или навигационного прием- ника, работающего с системой ГЛОНАСС, запишите данные координат и скорости ра- ботоспособного спутника и включите их в отчет по лабораторной работе. 4. Откройте файл ОRBITA_1.m из папки ОРБИТА_GLONASS_My ивыполните его. Ознакомьтесь с полученным графическим изображением орбиты спутника ГЛОНАСС и выполните задания 2, 3, 4. 5. Задание 2. Введите координаты и скорости спутника ГЛОНАСС из п. 3 в строки вход- ных данных файла ОRBITA_1.m, исполните файл. Проанализируйте полученное гра- фическое изображение орбиты спутника ГЛОНАСС. Результаты анализа внесите в от- чет. 6. Задание 3. Дополните файл ОRBITA_1.m процедурой вывода координат и скорости спутника на одно из значений текущего времени. Исполните файл и запишите полу- ченные значения координат и скоростей спутника в отчет. 7. Задание 4. Измените в п. 6 текущее время на 900 секунд граница интервала размноже- ния эфемерид), исполните файл занесите результат выполнения файла в отчет. Про- анализируйте и объясните разницу в результатах, полученных в п. п. 6 и 7. 4.2.3 Вопросы и задания для самоподготовки 1. Для каких целей требуется размножать координаты и скорости навигационных спут- ников? 2. Какие параметры являются начальными условиями при решении системы дифферен- циальных уравнении орбитального движения спутников ГЛОНАСС ? 3. В какой системе координат передаются данные о координатах и скорости в спутнико- вой системе ГЛОНАСС? 97 4. В какой системе координат решаются дифференциальные уравнения орбитального движения спутников ГЛОНАСС? 4.2.4 Файл ОRBITA_1.m %Имя m-файла:ORBITA_1.m %Программа иллюстрирует процедуру размножения эфемерид и орбиты спутника ГЛОНАСС %(демонстрация упрощенного варианта решения системы дифференциальных уравнений %движения спутника) %Программа выполняется совместно с функцией orbit_GL, использующей функцию MatLab ode45 %для решения дифференциальных уравнений методом Рунге-Кутта %Входные данные: %вектор координат x, y, z спутника XYZ (размерность-метр); %вектор скоростей спутника по осям x, y, z (размерность-м/с) VXYZ; %текущее время t= "начальное время" : "шаг" : "конечное время= "время в часах"*3600" %Выходные данные: %координаты спутника X, Y,Z (x, y, z) в абсолютной (относительной) системах координат; %скорости спутника Vx, Vy, Vz в абсолютной системе координат; %вектор текущего времени T; %вектор текущих координат и скоростей V %Расчет вектора входных параметров y omega = 0.7292115*10^(-4);%- скорость вращения Земли t=0:360:23*3600; S=-omega*3*3600;% угол %XYZ=[21840.10466;-9006.95351;-9696.59786];%координнаты спутника XYZ=[9795803.22265 ;-7174949.70703;22480344.23828 ];%координнаты спутника mS=[cos(S) -sin(S) 0;sin(S) cos(S) 0;0 0 1]; %матрица преобразования координат %VXYZ=[-1.19933288;0.58113958;-3.25131421];%скорости спутника VXYZ=[2773.857116;1295.602798;-814.5313262]; ys1=[mS*XYZ]';%вектор преобразованных координат ys2=[mS*VXYZ]'+omega*[-ys1(2) ys1(1) 0];% вектор преобразованных скоростей y=[ys1 ys2];%вектор начальных условий %Расчет орбиты спутника с помощью функции ode45 %[T,V] = ode45(@orbit_GL,[0:360:23*3600],[y],[]); [T,V] = ode45(@orbit_GL,[t],[y],[]); % Координаты и график орбиты спутника X=V(:,1); Y=V(:,2); Z=V(:,3) ; 98 subplot(2,1,1), plot3(X,Y,Z),grid on set(get(gcf,'CurrentAxes'),'FontSize',14,'FontName','TimesNewRoman'); %Vx=V(:,4); %Vy=V(:,5); %Vz=V(:,6) ; %subplot(1,3,2), plot3(Vx,Vy,Vz) % Координаты и график орбиты спутника в системе координат ПЗ90 S=omega*T; x= X.*cos(S)+Y.*sin(S); y =-X.*sin(S)+Y.*cos(S); z =Z; subplot(2,1,2), plot3(x,y,z),grid on set(get(gcf,'CurrentAxes'),'FontSize',14,'FontName','TimesNewRoman'); function [dy1 y1]= orbit_GL(t,y1) %Имя функции: orbit_GL %Функция записи системы дифференциальных уравнений для решения с помощью стандартной % программы MatLab dy1 = zeros(6,1); prom=398600.44*10^9/((y1(1)*y1(1)+y1(2)*y1(2)+y1(3)*y1(3))^1.5); dy1=[y1(4) y1(5) y1(6) [-y1(1) -y1(2) -y1(3)]*prom]'; 4.3 Орбитальное движение спутников ГЛОНАСС 4.3.1 Краткие сведения из теории Приведенный ниже комплекс программ для изучения орбитального движения спут- ников ГЛОНАСС составлен на основе формул для расчета орбит спутников по данным альманах ГЛОНАСС в полном соответствии с интерфейсным контрольным документом ГЛОНАСС [2]. На CD-диске комплекс размещен в папке 09_ORBITA_GL_NAVIOR. Структура комплекса изображена на рис. 4.7. Экспериментальные данные, использован- ные при изучении, поучены с помощью навигационного приемника «НАВИОР- 14», раз- работанного Государственным предприятием «Оризон- Навигация» (Украина). 99 Рис. 4.7. Структура комплекса программ орбитального движения спутников ГЛОНАСС Цель лабораторной работы: исследование орбитального движения спутников ГЛОНАСС по данным альманаха ORBITA_GLONASS.m Входные данные 'In_dat\GLN_8_11.alm' Функция «Чтение альманаха» Read_GL_Alm Функция преобразования дня при- вязки данных альманаха Gln_data_from_NA Функция JD_data Функция LLH_to_ECEF Функция WGS84_to_PZ90 Функция ECEF_to_LLH Функция GLN_satfind Функция ris_vis_sat Функция JD_epohi Функция init_data Функция JD_from_epohi Функция s0_Nut Функция llh_to_eci Функция init_satpos_gln Функция eci_to_ecef Функция top_coord Функция rewrite_satpos Функция gln_a_1 Функция satpos_eci_in_metr Функция init_satpos_gln Функция semi_axis_1 Функция kepler Функция utc_nut Функция utc_nut_fi_ep Функция koef Выходные данные: графики /данные 4.3.2 Лабораторная работа 4. 4 «Орбитальное движение спутников ГЛОНАСС» Рекомендуется следующий порядок выполнения лабораторной работы. 100 1. Создайте папку ORBITA_GL_NAVIOR_My и скопируйте в ее из папки ORBITA_GL_NAVIOR m-файлыи функции 2. Изучите функции и файл, используя рекомендованную литературу, комментарии к программам и блок –схему, изображенную на рис. 4.7. 3. Откройте m-файлОRBITA_1.m и выполните задания 1- 3, руководствуясь коммента- риями, приведенными в файле. 4. Задание 1. Постройте график орбит спутников ГЛОНАСС системе ECEF, проанализи- руйте их и результаты занесите в отчет. 5. Задание 2. Постройте график орбит спутников ГЛОНАСС системе ECI, проанализи- руйте их и результаты занесите в отчет. 6. Задание 3. Постройте графики времени наблюдения спутников ГЛОНАСС , проанали- зируйте их и результаты занесите в отчет. 4.3.3 Вопросы и задания для самоподготовки 1. Объясните, какой смысл вкладывается в содержание составляющих альманаха ГЛО- НАСС: поправка к шкале времени ГЛОНАСС относительно UTC(SU), номер четырех- летнего периода, поправка на расхождение системных шкал времени GPS и ГЛО- НАСС, календарный номер суток внутри четырехлетнего периода, номер спутника, номер несущей частоты, долгота восходящего узла орбиты спутника , время прохож- дения восходящего узла орбиты спутника, поправка к среднему значению наклонения орбиты спутника, поправка к среднему значению драконического периода обращения спутника, скорость изменения драконического периода обращения спутника, эксцен- триситет орбиты спутника, аргумент перигея орбиты спутника, признак состояния спутника. 2. Какая размерность данных передаваемых со спутника ГЛОНАСС в альманахе? 3. В какой системе координат передаются данные со спутника ГЛОНАСС. 4.3.4 Функции и файлы из папки ORBITA_GL_NAVIOR Файл ORBITA_GLONASS.m clear all %Имя файла:ORBITA_GLONASS.m %Программа рассчитывает орбиты спутников ГЛОНАСС по данным альманаха приемника НАВИОР 14 % разработанного ГП ОРИЗОН-НАВИГАЦИЯ в строгом соответствии с интерфейсным контрольным %документом ГЛОНАСС % Входные данные: 101 %Dat - альманах ГЛОНАСС находится в папке In_dat, например, Dat = 'In_dat\GLN_8_11.alm'; %Data_observ.year - год ; %Data_observ.mon - месяц ; %Data_observ.day - день; %current_loc_wgs - координаты точки, из которой проводится наблюдение спутников (pos. reciver (WGS-84): % current_loc_wgs.lat - широта, радиан ( DEG_TO_RAD * (50. + 29.0 / 60.0 + 36.78 / 3600.0),например, % 50 градусов; 29 минут; 36.78 секунд; % current_loc_wgs.lon - долгота, радиан ( DEG_TO_RAD * (30. + 27.0 / 60.0 + 50.5 / 3600.0),например, %30 градусов; 27 минут; 50.5 секунд; %current_loc_wgs.h - вісота, метр, например, 120.9;все в WGS 84 %step_t - шаг в секундах (600= 10 минут); %L =24*1 * 3600- расчетный интервал, например, 24 часа, 1 день, 3600 секунд %kol - количество спутников, для которых строятся графики орбитального движения %nom_ns(1:kol) - номера спутников, для которых строятся графики орбитального движения (записы- ваются % в квадратных скобках, через пробел, число номеров должно совпадать с количеством спутников % %Выходные данные: %x(j,i)= satpos_gln_ecef(i).x - координата спутника x c номером i на момент времени j в системе ECEF; % y(j,i)= satpos_gln_ecef(i).y- координата спутника у c номером i на момент времени j в системе ECEF; %z(j,i)= satpos_gln_ecef(i).z- координата спутника z c номером i на момент времени j в системе ECEF; %x1(j,i)= satpos_gln_a(i).x- координата спутника x c номером i на момент времени j в системе ECI; %y1(j,i)= satpos_gln_a(i).y- координата спутника y c номером i на момент времени j в системе ECI; %z1(j,i)= satpos_gln_a(i).z- координата спутника z c номером i на момент времени j в системе ECI; %x2(j,i)= satpos_gln_ecef(i).vx - cкорость спутника c номером i на момент времени j по оси x; %y2(j,i)= satpos_gln_ecef(i).vy- cкорость спутника c номером i на момент времени j по оси y; %z2(j,i)= satpos_gln_ecef(i).vz- cкорость спутника c номером i на момент времени j по оси z; %x3(j,i)= satvis_gln(i).el*180/pi - угол места спутника c номером i на момент времени j; %y3(j,i)= satvis_gln(i).az- угол азимута спутника c номером i на момент времени j; %z3(j,i)= satvis_gln(i).r - дальность до спутника c номером i на момент времени j; %plot3(x(:,prn),y(:,prn),z(:,prn),S,'LineWidth',0.5)- график орбит спутников ГЛОНАСС системе ECEF; %plot3(x1(:,prn),y1(:,prn),z1(:,prn),S, 'LineWidth',1)- график орбит спутников ГЛОНАСС системе ECI; %ris_vis_sat- график видимости спутников %Примечание: места ввода данных отмечены строками %%%%%%%%%%%%%%%%%%%%% %Ввод входных данных Dat = 'In_dat\GLN_8_11.alm';%%%%%%%%%%%%%%%%%%%%% N=6378136;% радиус Земли (используется, как нормирующий коэффициент map(N);%функция выводит на графики Землю 102 %Задание цветов для графики j_color = 0; color6(1:11) = ['B' 'r' 'D' 'c' 'g' 'k' 'm' '.' 's' 'H' '+']; %color6(1:16) = [':' 'k' '.' 'r' 'g' 'r' 'c' 'm' 'r' ':' 'g' ':' 'b' ':' 'k' 'h']; [alm,max_kol] = Read_GL_Alm(Dat); % Чтение альманаха nom = 1; day_from_leap = alm(nom).Na; while (alm(nom).Health > 0 ) day_from_leap = alm(nom).Na; % номер дня от ближайшего предшествующего високосного года nom = nom + 1; end; leap_year = 2004;% високосный год %alm.Na -(сек) время привязки альманаха от начала предшествующего високосного года % Время привязки альманаха: timeUTC = Gln_data_from_NA(leap_year, day_from_leap); [JD_alm, day_year_alm] = JD_data(timeUTC); %JD_alm - номер юлианского дня привязки альманаха % Data_observ - Дата начала обсервации (расчета): Data_observ = timeUTC; %{ Data_observ.year = timeUTC.year ; Data_observ.mon = timeUTC.mon ; Data_observ.day = timeUTC.day + 2 ; %} Data_observ.year = 2006;%%%%%%%%%%%%%%%%%%%%% Data_observ.mon = 11; %%%%%%%%%%%%%%%%%%%%% Data_observ.day =9; %%%%%%%%%%%%%%%%%%%%% Data_observ.ti = 0 ; %%%%%%%%%%%%%%%%%%%%% [JD_observ, day_year_observ] = JD_data(Data_observ); % JD_observ - номер юлианского дня обсервации day_observ_alm = JD_observ - JD_alm; %вывод в командное окно даты, на которую выполняется расчет параметров (год, месяц, день) fprintf('year=%i mon = %i day = %i \n',Data_observ.year, Data_observ.mon, Data_observ.day); %Расчет количества и номеров здоровых спутников по данным альманаха: nom = 1; i = 0; k = 0; i = 0; 103 id = alm(nom).ID; while ( i < 24) id = alm(nom).ID; Health = alm(nom).Health; %fprintf('1: i=%i k=%i nom=%i id=%i Health=%i \n', i, k, nom, id, Health); if ( id > 0) Health = alm(nom).Health; if ( Health == 0) k = k + 1; nom_ns(k) = id; % fprintf('2: i=%i k=%i nom=%i id=%i Health=%i \n', i, k, nom, id, Health); nom = nom + 1; else nom = nom + 1; end; else nom = nom + 1; end; i = i + 1; end; % i kol = k; % fprintf('kol=%i \n', kol); % для вывода номеров здоровых спутников перед fprintf убрать "%" nom_ns % - номера здоровых навигационных спутников, %количество спутников, для которых строятся графики орбитального движения kol = 1;%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nom_ns(1:kol)=[24];%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %nom_ns(1:kol)=[1 2 3 4 7 8 19 21 22 23 24 ];%номера спутников (количество номеров должно %равняться количеству спутников KOL_GLN = 24; DEG_TO_RAD = 0.017453292519943; % (PI / 180.00 ) RAD_TO_DEG = 57.295779513082; % (180.00 / PI ) A_WGS84_M = 6378137.0 ; % WGS-84 ellipsoid parameters B_WGS84_M = 6356752.314; % WGS-84 ellipsoid parameters A_PZ90_M = 6378136.0; %6 378 136 м - Equatorial radius of the Earth - ,большая полуось эллипсоида B_PZ90_M = 6356751.36174571344; %AP_LAND (m) Polar radius of the Earth %A_PZ90_KM = 6378.136; %(Km) Equatorial radius of the Earth B_PZ90_KM = 6356.75136174571344; % AP_LAND (Km) Polar radius of the Earth FACTOR_PZ90 = 1.0/298.257839303;% Коэффициент сжатия эллипсоида C_LIGHT_M = 2.99792458E8; % m/sec Speed of light RAD_IN_SEC = 7.2722052166430e-5 ;%(PI/43200.0) Number radian in second of time 2*PI/(24*3600)=PI/43200; 104 dt_lsf = 14; % Вводим координаты приемника: % current_loc_wgs - Pos. reciver (WGS-84): current_loc_wgs.lat = DEG_TO_RAD * (50. + 29.0 / 60.0 + 36.78 / 3600.0);%%%%%%%%%%%%%%%% current_loc_wgs.lon = DEG_TO_RAD * (30. + 27.0 / 60.0 + 50.5 / 3600.0);%%%%%%%%%%%%%%%% current_loc_wgs.h = 120.9; % метр, WGS 84 %%%%%%%%%%%%%%%% %current_loc_wgs.hae = 0.1229; % Km PZ90 %Преобразование координат rec_pos_xyz_wgs = LLH_to_ECEF( A_WGS84_M, B_WGS84_M, current_loc_wgs); rec_pos_xyz_pz90_m = WGS84_to_PZ90(rec_pos_xyz_wgs); current_loc_pz90 = ECEF_to_LLH(A_PZ90_M, B_PZ90_M, rec_pos_xyz_pz90_m); step_t = 600;%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %L = 12*24; L =24*1 * 3600;%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %ti_start = 0; ti_start = Data_observ.ti + day_observ_alm * 86400; ti_end = ti_start + L; j = 0; for ( ti = ti_start : step_t : ti_end ) j = j + 1; %ПЕРЕСЧЕТ ВРЕМЕНИ из секунд в часы, минуты и секунды: t(j) = ti; A(j) = mod(t(j), 86400); hour(j) = floor(A(j) / 3600); m(j) = floor((A(j) - hour(j)*3600) / 60); sek(j) = A(j) - hour(j) * 3600 - m(j) * 60; for (n=1:24) % обнуление массивов в текущий момент времени для всех спутников x(j, n) = 0; y(j, n) = 0; z(j, n)= 0; ris_vid(n, j) = 0; end; %timeUTC.ti = mod( ti, 86400); % текущее время обсервации от начала суток timeUTC.ti = ti; 105 [kol_gln_a, satvis_gln, satpos_gln_a, satpos_gln_ecef] = GLN_satfind( A_PZ90_M, B_PZ90_M, timeUTC, current_loc_pz90, alm); for (i=1: KOL_GLN) if ( satvis_gln(i).el >= 0) el = RAD_TO_DEG * satvis_gln(i).el; az = RAD_TO_DEG * satvis_gln(i).az; AZ(j,i) = az;%уголы азимута спутников в градусах EL(j,i) = el;%уголы видимости спутников в градусах ris_vid(i, j) = i; % fprintf('i=%2i el=%6.2f az=%7.2f x=%12.2f y=%12.2f z=%12.2f \n', i, el, az, satpos_gln_a(i).x, sat- pos_gln_a(i).y, satpos_gln_a(i).z); else AZ(j,i) = 0.0; EL(j,i) = 0.0; end; end; for (i=1:KOL_GLN) prn = alm(i).ID; health = alm(i).Health; if ( (prn > 0) & (health == 0)) %Координаты спутников в системе ECEF. Для визуализации добавить в вывод графиков и % снять блокировку. x(j,i)= satpos_gln_ecef(i).x; y(j,i)= satpos_gln_ecef(i).y; z(j,i)= satpos_gln_ecef(i).z; %Координаты спутников в системе ECI. Для визуализации добавить в вывод графиков. x1(j,i)= satpos_gln_a(i).x; y1(j,i)= satpos_gln_a(i).y; z1(j,i)= satpos_gln_a(i).z; %Скорости спутников. Для визуализации добавить в вывод графиков. x2(j,i)= satpos_gln_ecef(i).vx; y2(j,i)= satpos_gln_ecef(i).vy; z2(j,i)= satpos_gln_ecef(i).vz; %Углы видимости и дальности до спутников. Для визуализации добавить в вывод графиков. x3(j,i)= satvis_gln(i).el*180/pi; y3(j,i)= satvis_gln(i).az; z3(j,i)= satvis_gln(i).r; end; end; end; % ti ( j ) 106 max_n = 24; %для блокировки вывода графика "Время наблюдения спутников ГЛОНАСС" перед функцией %ris_vis_sat в следующей строке поставить % ris_vis_sat(max_n, j, ti_start, step_t, ris_vid); for (i=1:kol) j_color = j_color + 1; if (j_color > 14 ) j_color = 1; end S = color6(j_color); prn = nom_ns(i); hold on h_F1 = gca; %{ %ГРАФИКИ ПАРАМЕТРОВ ОРБИТАЛЬНОГО ДВИЖЕНИЯ СПУТНИКОВ ГЛОНАСС %График орбит спутников ГЛОНАСС системе ECEF plot3(x(:,prn),y(:,prn),z(:,prn),S,'LineWidth',0.5); axis([ -2.552*10^(7) 2.552*10^(7) -2.552*10^(7) 2.552*10^(7) -2.552*10^(7) 2.552*10^(7)]); set(get(gcf,'CurrentAxes'),'FontSize',14,'FontName','TimesNewRoman'); set(h_F1,'Position',[0.1 0.1 0.85 0.9]) ; xlabel('Координата X') ylabel('Координата Y'), zlabel ('Координата Z'),grid on str1 = num2str( prn); text(x(j,prn), y(j,prn),z(j,prn),str1,'FontSize',14,'FontName','TimesNewRoman','HorizontalAlignment','center' ); %} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %{ %График орбит спутников ГЛОНАСС системе ECI plot3(x1(:,prn),y1(:,prn),z1(:,prn),S, 'LineWidth',1); axis([ -2.552*10^(7) 2.552*10^(7) -2.552*10^(7) 2.552*10^(7) -2.552*10^(7) 2.552*10^(7)]) set(get(gcf,'CurrentAxes'),'FontSize',14,'FontName','TimesNewRoman'); set(h_F1,'Position',[0.1 0.1 0.85 0.9]) ; xlabel('Координата X') ylabel('Координата Y'), zlabel('Координата Z'),grid on str1 = num2str( prn); 107 |