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

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

Голосов: 113

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

Приведенный ниже текст получен путем автоматического извлечения из оригинального PDF-документа и предназначен для предварительного просмотра.
Изображения (картинки, формулы, графики) отсутствуют.
    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



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