Единое окно доступа к образовательным ресурсам

Спутниковые системы навигации: Лабораторный практикум на компьютере

Голосов: 113

Приводится лабораторный практикум по вопросам спутниковых радионавигационных систем в виде открытых программных комплексов в среде MatLab. Лабораторный практикум в виде 16 работ с заданиями, примерами и контрольными вопросами охватывает направления формирования сигналов спутниковых навигационных систем GPS и ГЛОНАСС, преобразования координат, моделирование орбит навигационных спутников GPS, ГЛОНАСС, GALILEO, декодирования и расшифровка данных навигационных спутников, решение навигационной задачи расчета позиции приемника пользователя. Для студентов, аспирантов и преподавателей дисциплин, изучающих системы и компоненты спутниковых средств навигации.

Приведенный ниже текст получен путем автоматического извлечения из оригинального PDF-документа и предназначен для предварительного просмотра.
Изображения (картинки, формулы, графики) отсутствуют.
    3. Запустите MatLab [7, 8].
4. Обратитесь к папке ORBITA_GPSv1_My и откройте ее.
5. Последовательно откройте функции и m- файлы из папки ORBITA_GPSv1:
   Orbita_GPS.m, Yuma_GPS_Alm1.m, Tim.m, LLHECEF.m, ECEFLLH.m, PR_Tim.m,
   PR_Yuma_GPS_Alm1.m, almanac_yuma_week0371_589824.txt. Изучите тексты про-
   грамм и комментарии в каждой программе. Выполните задания 2 и 3.
6. Задание 2. Выполните m- файл PR_Tim.m и изучите в командном окне полученный
   результат. Сформируйте на основе m- файла PR_Tim.m файл с именем PR_Tim_My.m.
   Введите в сформированный файл данные, соответствующие дате выполнния работы,
   сохраните файл и выполните файл. Результат из командного окна запишите в отчет.
7. Задание 3. . Выполните m- файл PR_Yuma_GPS_Alm1.m и изучите в командном окне
   полученный результат. Сопоставьте данные из командного окна с данными альманаха
   almanac_yuma_week0371_589824.txt.     Сформируйте      на    основе     m-   файла
   PR_Yuma_GPS_Alm1.m файл с именем PR_Yuma_GPS_Alm1_My.m. Введите в сфор-
   мированный файл имя записанного альманаха, сохраните файл и выполните файл,
   сравните результат из командного окна с данными альманаха и результат запишите в
   отчет.
8. Очистите командное окно из папки ORBITA_GPSv1_My, откройте m файл с именем
   Orbita_GPS.m, выполните m-файл, ознакомьтесь с графическим результатом и выпол-
   ните задание 4.
9. Задание 4. Присвойте переменной Dat в файле Orbita_GPS.m имя альманаха из п.1 и
   введите данные о начале отсчета как указано в комментариях к файлу. Установите
   скорость вращения Земли равную нулю. Убедитесь, что все графики, кроме графика 1
   закомментированы. Исполните файл. Результат выполнения график 1 запишите на
   дискету, как приложение к отчету. В отчет внесите описание графика 1.
10. Задание 5. Установите скорость вращения Земли, равную 7.2921151467e-005. Выбери-
   те номер спутник, с которым будете работать. Рекомендация: для генерации вариантов
   для каждого студента можно задавать номер спутника равный, например, номеру дня
   рождения, поскольку число спутников GPS, как правило, близко к 30. Разблокируйте
   график 2, исполните файл. Дайте объяснения зависимостям улов азимута и видимости.
   График запишите на сменный диск, а интерпретацию результатов в отчет.
11. Задание 5. Разблокируйте график 3, исполните файл. Дайте объяснения зависимостям
   доплеровской частоты, широты от долготы. График запишите на сменный диск, а ин-
   терпретацию результатов в отчет.



                                          71


12. Задание 6. Разблокируй график 4, исполните файл. Объясните поведение проекции
   орбиты на плоскость XY и изменение дальности до спутник от времени. График запи-
   шите на сменный диск, а интерпретацию результатов в отчет.
13. Задание 7. Выберите группу спутников, принадлежащих одной из 6 орбит и повторите
   расчетные процедуры заданий 4- 6. График запишите на сменный диск, а интерпрета-
   цию результатов в отчет.
14. Оформите отчет. Графики результатов моделирования можно предъявить на сменном
   диске, как приложение к отчету.

      4.1.3 Вопросы и задания для самоподготовки

1. Объясните, какой смысл вкладывается в содержание составляющих альманаха YUMA:
   ID, Health, Eccentricity, Time of Applicability(s), Orbital Inclination(rad), Rate of Right As-
   cen(r/s), SQRT(A) (m 1/2), Right Ascen at Week(rad), Argument of Perigee(rad), Mean
   Anom(rad), Af0(s) Af1(s/s), week.
2. Объясните зависимость изменения доплеровской частоты при движении спутника по
   орбите.
3. Объясните зависимость изменения дальности до спутника от времени для неподвиж-
   ного наблюдателя.
4. Для каких целей используются данные альманаха в спутниковых навигационных при-
   емниках.
5. Запишите уравнение для расчета дальностей до спутников и найдите это уравнение в
   текстах программ.
6. С помощью какого фрагмента программного комплекса рассчитываются углы видимо-
   сти и азимута спутников?
7. Какой содержательный смысл имеет таблица в командном окне, появляющаяся после
   выполнения файла Orbita_GPS.m.

      4.1.4 Функции и файлы из папки ORBITA_GPSv1


     Функция ECEFLLH
     function [Rx,Ry,Rz] = ECEFLLH(lon, lat,hr)
     %Имя функции: ECEFLLH.m
     %Функция выполняет преобразование координат
     %Входные данные:lon-долгота,lat-широта,h-высота;a,b-большая
     %и малая полуоси эллипсоида
     %Выходные данные:Rx,Ry,Rz- координаты в ECEF
     %Для WGS-84


                                                  72


a=6378137.0;
b=6356752.314;
 n=a*a/sqrt(a*a*cos(lat)*cos(lat)+b*b*sin(lat)*sin(lat));
 Rx=(n+hr)*cos(lat)*cos(lon);
 Ry=(n+hr)*cos(lat)*sin(lon);
 Rz=(b*b/(a*a)*n+hr)*sin(lat);
Функция LLHECEF
function [lons,lats,hrs] = LLHECEF(Xk,Yk,Zk)
%Имя функции: LLHECEF.m
%Функция выполняет преобразование координат.
%Входные данные:Rx,Ry,Rz- координаты в ECEF
%Выходные данные:lon-долгота,lat-широта,h-высота
%a,b-большая и малая полуоси эллипсоида
a=6378137.0;
b=6356752.314;
  xy = sqrt(Xk*Xk + Yk*Yk);
  thet = atan(Zk*a/(xy*b));
  esq = 1.0-b*b/(a*a);
 epsq = a*a/(b*b)-1.0;
  lats = atan((Zk+epsq*b*(sin(thet)^3))/(xy-esq*a*(cos(thet)^3)));
  lons = atan2(Yk,Xk);%!
if lons < 0
  lons = 2*pi + lons;
end ;
 n = a*a/sqrt(a*a*cos(lats)*cos(lats) + b*b*sin(lats)*sin(lats));
  hrs = xy/cos(lats)-n;
end
Функция Tim
function [week,modeweek,d,dweek,weeks]=Tim(d2,h,min,s)
%Имя функции:Tim.m
%Функция Tim.m работает совместно со стандартной функцией MatLab DAYSDIF.m и рассчитывает:
%week- текущую GPS-неделю, modeweek- модифицированную GPS неделю, d- количество дней,
%dweek- день недели ,weeks- время GPS
%Входные данные:d2='10/23/2007' - 'номер месяца/номер дня месяца/номер года', h=23.0 -часы,
%;min=59.0- минуты, s=59.0- секунды на которые рассчитываются выходные данные
%d2='10/35/2003';h=23.0;min=59.0;s=59.0;
week = floor(DAYSDIF('1/6/1980',d2,3)/7);% текущая GPS-неделя
modeweek=week-1024;% модифицированная GOS-неделя
d = DAYSDIF('1/6/1980',d2,3);%количество дней
dweek=fix(d-week*7);%номер дня недели (нулевой день-воскресенье)
weeks=(dweek)*24*60*60+h*60*60+min*60+s;% время GPS в неделе (секунды)
Файл PR_Tim.m
                                                  73


       %Пример применения функции Tim.m
       %Имя m- файла: PR_Tim.m
       d2= '10/13/2006';
       h=22.0;
       min=40.0;
       s=11.0;
       [week,modeweek,d,dweek,weeks]=Tim(d2,h,min,s);
       [week,modeweek,d,dweek,weeks] %=1396             372      9777     5      513611-результат расчета;
       %1396-неделя GPS, отсчитываемая с ночи с 5 на 6 января 1980 года, 372=1396-1024- модифициро-
ванная
       % неделя GPS, 9777- количество дней прошедших с 6 января 1980 года, 5-пятый день недели (пятни-
ца),
       %считая с понедельника, 513611-количество секунд от начала текущей недели.
       d = DAYSDIF('1/6/1980',d2,3);%функция MatLab
       Функция Yuma_GPS_Alm1
       function [alm,max_kol] = Yuma_GPS_Alm1(Dat)
       %Имя функции:Yuma_GPS_Alm1.m
       %Функция читает данные альманаха, записанные в формате YUMA
       %Входные данные записываются в переменную Dat, например,
       %Dat='Имя файла альманаха в формате YUMA'
       %Выходные данные:1. Численные значения альманаха спутников GPS, представляемые
       %в виде структуры в переменной alm =[%alm(ID).ID(1); alm(ID).Health(2); alm(ID).e(3);
       %alm(ID).TOA(4); alm(ID).deltai(5);%alm(ID).OMEGADOT(6); alm(ID).A05(7); alm(ID).omega0(8);
       %alm(ID).omega(9);%alm(ID).M0(10); alm(ID).Af0(11); alm(ID).Af1(12); alm(ID).Week(13)], где
       % цифра в скобках обозначает порядковый номер параметра альманаха в формате YUMA .
       %Для       чтения      альманаха   в   m-файл    фикция    записывается     в   виде   [alm,max_kol]   =
Yuma_GPS_Alm1(Dat).
       %2. Количество спутников GPS записывается в переменную max_kol.
       for i=1:32% цикл
         alm(i).ID = 0;% обнуление массива
       alm(i).Health=63;% обнуление массива
       end;


       fid =fopen(Dat,'rt');% открыть файл для чтения
       %чтение данных из файла
       max_kol = 0;
       while not(feof(fid))


          s1=fscanf(fid,'%s',6);
          if not(feof(fid))
              lenstr = length(s1);



                                                        74


      while (fscanf(fid,'%s',1) == '*')
      end
      str1 = fscanf(fid,'%s',1);
      lenstr = length(str1);
      n_sv = sscanf(str1,'%d');
      strID = str1(1:lenstr);
      ID = sscanf(strID,'%d');
      alm(ID).ID = ID;


      t_2=fscanf(fid,'%s',1);
      alm(ID).Health=fscanf(fid,'%d',1);


      t_3=fscanf(fid,'%s',1);
      alm(ID).e = fscanf(fid,'%g',1);


      t_4=fscanf(fid,'%s',3);
      alm(ID).TOA =fscanf(fid,'%g',1);
      t_5=fscanf(fid,'%s',2);
      alm(ID).deltai=fscanf(fid,'%g',1);%i0
      t_6=fscanf(fid,'%s',4);
      alm(ID).OMEGADOT=fscanf(fid,'%g',1);
      while not(fscanf(fid,'%c',1) == ':')
      end
      alm(ID).A05=fscanf(fid,'%g',1);
      t_8=fscanf(fid,'%s',4);
      alm(ID).omega0 =fscanf(fid,'%g',1);
      t_9=fscanf(fid,'%s',3);
      alm(ID).omega=fscanf(fid,'%g',1);
      t_10=fscanf(fid,'%s',2);
      alm(ID).M0=fscanf(fid,'%g',1);
      t_11=fscanf(fid,'%s',1);
      alm(ID).Af0=fscanf(fid,'%g',1);
      t_12=fscanf(fid,'%s',1);
      alm(ID).Af1=fscanf(fid,'%g',1);
      t_13=fscanf(fid,'%s',1);
      alm(ID).Week=fscanf(fid,'%g',1);
      max_kol = max_kol+1 ;
  end
end
fclose(fid)
Файл Orbita_GPS.m
clear all

                                              75


     %Имя m-файла:Orbita_GPS.m
     %Программа рассчитывает орбиты навигационных спутников GPS
     %Входные данные:
     %файл альманаха в формате Yuma,имя файла альманаха присваивается
     %переменной "Dat",например,Dat = 'имя файла альманаха';
     %данные о начале отсчета "d2",d2='месяц/день/год';h=час;min=минута;s=секунда;
     %координаты позиции приемника –lat (широта в радианах),lon (долгота в радианах,
     %hr (высота в метрах);
     %шаг, с каким будут рассчитываться параметры орбит (step,секунды);
     %количество точек (L), в которых будут рассчитываться параметры орбит
     %L=12*3600/step,L читается так: количество часов (например,12)
     %число секунд в часе (3600) деленное на шаг (step)
     %В программе применяются функции: Yuma_GPS_Alm1.m- считывание данных альманаха,
     %заданного в формате YUMA; ECEFLLH.m, LLHECEF.m - преобразование координат;Tim.m- расчет
времени;
     %Постоянные:
     %скорость вращения Земли
     OMEGAeDOT=7.2921151467e-005;
     %или
     %OMEGAeDOT=0;
     mu=398600500000000;
     F_CONST                    = 4.442807633E-10;
     %Задание цветов для графики
     j_color = 0;
     color6(1:16) = [':' 'k' '.' 'r' 'g' 'r' 'c' 'm' 'r' ':' 'g' ':' 'b' ':' 'k' 'h'];
     %Входные данные
     Dat = 'almanac_yuma_week0371_589824.txt';
     d2='10/06/2006';h=13.0;min=8.0;s=55.0;
     lat = 0.88032730015257; %50 град; 26 мин.; 20.54 с
     lon = 0.53109641675259;%30 град; 25 мин.; 46.4995 с
     hr=184;%высота в метрах
     step=300;
     L=(10*3600)/step;


     %Чтение альманаха


     [alm,max_kol] = Yuma_GPS_Alm1(Dat);
     kol = 0;
     for i = 1 : max_kol
     id=alm(i).ID;
     if id > 0
     kol = kol + 1;


                                                                          76


nom_ns(kol) = id;
end
nom_ns;%номер навигационного спутника
end


%Преобразование координат


[Rx,Ry,Rz] = ECEFLLH(lon, lat,hr);
%Rx=0;Ry=0;,Rz=0;%центр масс Земли
%Выбор спутников:
%для выбора спутников вводится параметры kol-количеество спутников для
%исследования и номера спутников, например, kol =4; nom_ns(1:kol) = [3 6 7 31],
%такая запись обозначает, что исследуются (рассчитываются орбиты 4 спутников
%с номерами 3,6,7,31; количество номеров спутников должно совпадать с kol
%Варианты (можно любые другие)
%kol =9
%nom_ns(1:kol) = [1 3 4 5 6 7 8 9 10];
%kol =5
%nom_ns(1:kol) = [1 13 14 26 29];        %1 спутники орбиты 1
%kol =5
%nom_ns(1:kol) = [2 5 22 28 30];         %2 спутники орбиты 2
%kol =4
%nom_ns(1:kol) = [3 6 7 31];         %3 спутники орбиты 3
%kol =5
%nom_ns(1:kol) = [4 11 15 17 24 ];       %4 спутники орбиты 4
%kol =4
%nom_ns(1:kol) = [8 9 25 27];        %5 спутники орбиты 5
%kol =5
%nom_ns(1:kol) = [10 18 20 21 23];        %6 спутники орбиты 6
%kol =29
%nom_ns(1:kol) = [1 3 4 5 6 7 8 9 10 11 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31] ;
%kol =14;
%nom_ns(1:kol) = [1 3 4 5 6 7 8 9 10 11 13 14 15 16];% 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31] ;


kol =2;
nom_ns(1:kol)=[1 3 ];
str1 = num2str(nom_ns(1:kol));
for k = 1 : kol
 i = nom_ns(k);
%Начало отсчета текущего времени
[week,modeweek,d,dweek,weeks]=Tim(d2,h,min,s);
%Расчет орбит спутников по формулам ( )


                                                 77


for j = 1:L % 0:L
   t( j )=weeks+step*j; %-step;
  %t1(j) = t(j)/60; %изменение дискретности текущего времени
  %d_wn = (week - alm(i).Week);%если в альманахе учтены 1024
 %d_wn = 0;
  d_wn=(modeweek-alm(i).Week);%если в альманахе не учтено 1024
  tk = t(j) + d_wn * 604800 - alm(i).TOA;
  d_wn = abs(modeweek - alm(i).Week);
 dd = 302400.0 + d_wn * 604800;
while (abs(tk) > dd)
  if tk > dd
  tk = tk - 604800;
  else
     if tk < -dd
         tk = tk + 604800;
            end
  end % if
end % while


%Справочник по альманаху- цифра в скобках обозначает порядковый номер
%параметра альманаха в формате YUMA
%alm(ID).ID(1); alm(ID).Health(2); alm(ID).e(3); alm(ID).TOA(4); alm(ID).deltai(5);
%alm(ID).OMEGADOT(6); alm(ID).A05(7); alm(ID).omega0(8); alm(ID).omega(9);
%alm(ID).M0(10); alm(ID).Af0(11); alm(ID).Af1(12); alm(ID).Week(13);


  n0=sqrt((mu)/(alm(i).A05^6));
  n=n0;
  Mk = alm(i).M0+n*tk;
  e=alm(i).e;


  %Решение уравнения Кеплера
 eps = 1.0E-15;
 y = e * sin(Mk);
 x1 = Mk;
 x = y;
 for k = 0 : 15
    x2 = x1;
    x1 = x;
    y1 = y;
    y = Mk - (x - e * sin(x));
    if (abs(y - y1) < eps)
         break


                                               78


         end
   x = (x2 * y - x * y1) / (y - y1);
  end %k
  Ek = x;
   deltr = F_CONST * alm(i).e * alm(i).A05 * sin(Ek);
 dt1 = alm(i).Af0 + alm(i).Af1 * tk + deltr;
 tk = tk - dt1;
nuk =atan2(sqrt(1-alm(i).e^2)*sin(Ek),(cos(Ek)-alm(i).e));
Ek = acos((alm(i).e+cos(nuk))/(1+alm(i).e*cos(nuk)));
Fk =nuk + alm(i).omega;
uk =Fk;
ik=alm(i).deltai;
rk =(alm(i).A05^2)*(1.0-alm(i).e*cos(Ek));
xkk =rk*cos(uk);
ykk =rk*sin(uk);
OMEGAk =alm(i).omega0+(alm(i).OMEGADOT-OMEGAeDOT)*tk-OMEGAeDOT*alm(i).TOA;
%Координаты спутников
Xk(j) = xkk *cos(OMEGAk)-ykk*cos(ik)*sin(OMEGAk);
Yk(j) = xkk*sin(OMEGAk)+ykk*cos(ik)*cos(OMEGAk);
Zk(j) = ykk*sin(ik);
%Дальности до спутников
PR(j) = sqrt((Xk(j) - Rx)^2 + (Yk(j) - Ry)^2 + (Zk(j) - Rz)^2);


%Перевод в географическую систему если требуется
%[lons,lats,hrs] = LLHECEF(Xk,Yk,Zk);
%(Llon(j),Llat(j),Hhr(j)) = [lons,lats,hrs];


xls = Xk(j) - Rx;
yls = Yk(j) - Ry;
zls = Zk(j) - Rz;
range1 = sqrt(xls*xls+yls*yls+zls*zls);
if j>1
doppler(j-1) = (range1 - range2) * 5.2514 / step;%расчет доплеровской частоты
end
range2 = range1;
P = sqrt(Rx * Rx + Ry * Ry + Rz * Rz);
tdot = ( Rx*xls+Ry*yls+Rz*zls)/range1/P;
xll = xls/range1;
yll = yls/range1;
zll = zls/range1;
%Расчет угла видимости
if tdot >= 1.00


                                                 79


b = 0.0;
elseif tdot <= -1.00
b = pi;
else
b = acos( tdot);
end
satang = pi/2.0 - b;
TT(j) =satang;


%Расчет угла азимута
xn =-cos(lon)*sin(lat);
yn =-sin(lon)*sin(lat);
zn = cos(lat);
xe =-sin(lon);
ye = cos(lon);
xaz = xe*xll + ye*yll;
yaz    = xn*xll + yn*yll + zn*zll;
if (xaz == 0) or (yaz == 0)
az(j)= 0;
else
az(j) = atan2(xaz,yaz);
end
if az(j) < 0
az(j) = az(j) + pi*2;
end
end % j
for j = 1:L
[Llon(j),Llat(j),Hhr(j)] = LLHECEF(Xk(j),Yk(j),Zk(j));%преобразование координат
if j > 1
if abs(Llon(j)-Llon(j-1)) > pi
Llon(j) = Llon(j) + 2*pi;
       end
  end
end


j_color = j_color + 1;
if (j_color > 14 )
j_color = 1;
end
%F_ont=get(gcf,'CurrentAxes'),'FontSize',16,'FontName','TimesNewRoman';%формат текста на графиках
S = color6(j_color);
%Графика


                                            80



    
Яндекс цитирования Яндекс.Метрика