Спутниковые системы навигации. Учебное пособие. Учебное пособие ( Лабораторный практикум на компьютере ) Київ 2008 1
 Скачать 4.49 Mb.
|
|
Функция CRC32Value function[ulCRC] = CRC32Value(i) %Имя функции: CRC32Value %Функция для контроля декодированных данных CRC32_POLYNOMIAL = hex2dec('EDB88320'); ulCRC = i; for (j = 1:8) a = bitand(ulCRC, 1); if a > 0 ulCRC = bitxor(bitshift(ulCRC, -1), CRC32_POLYNOMIAL); else ulCRC = bitshift(ulCRC, -1); end; end; 171 Функция CalculateBlockCRC32 function [ulCRC] = CalculateBlockCRC32(fid, size1) %Имя функции:CalculateBlockCRC32 %Функция для проверки CRC-кода ulCRC = 0; int_ulcrc = uint(32); FF0 = hex2dec('00ffffff'); ff = hex2dec('ff'); ulCount = size1; while ( ulCount > 0 ) ulTemp1 = bitand( bitshift(ulCRC, -8 ), FF0); ucBuffer = fread(fid,1,'uint8'); int_ulcrc = ulCRC; xor1 = bitxor(int_ulcrc, ucBuffer ); param1 = bitand(xor1, ff); ulTemp2 = CRC32Value(param1 ); ulCRC = bitxor(ulTemp1, ulTemp2); ulCount = (ulCount - 1); end Функция e_norm function str_norm = e_norm(e_dat, mantissa); %Имя функции:e_norm %Функция для форматирования записи данных format_e = num2str(mantissa); format_e = strcat('%0.',format_e,'E'); str_norm = num2str(e_dat*10,format_e); k=strfind(str_norm,'.'); if k>2 str_norm = strcat(str_norm(1),'0.',str_norm(2:k-1), str_norm(k+1:length(str_norm))); else str_norm = strcat(' 0.',str_norm(1:k-1), str_norm(k+1:length(str_norm))); end; Функция var_sign function [result] = var_sign(a, poz, koef) %Имя функции:var_sign %Вспомогательная функция для декодирования коэффициентов af0, af1 %af0_bin = dec2bin(a) sign_a = bitget(a,poz); 172 a_plus = bitset(a, poz, 0); %a_plus_bin = dec2bin(af0_plus) if sign_a==0 result = a_plus * koef; else a_minus = bitcmp(a_plus,(poz - 1)) + 1; %a_bin = dec2bin(a_minus) result = - a_minus * koef; end end 5.3 Модель движения навигационных спутников GPS и ГЛОНАСС 5.3.1 Краткие сведения из теории Развитие систем спутниковой радионавигации в настоящее время и в ближайшей перспективе ориентировано на применении нескольких систем. При этом в одном навига- ционном приемнике должны обрабатываться сигналы спутников систем, которые нахо- дятся в зоне видимости. В представленном комплексе программ средствами MatLab [7, 8] решена задача моделирования движения спутников GPS и ГЛОНАСС с применением дан- ных альманаха спутниковых систем. Программирование выполнено по алгоритмам пол- ностью соответствующим интерфейсному контрольному документу [9]. Программа рассчитывает углы азимута и места навигационных спутников GPS и ГЛОНАСС по данным совместного для обеих систем альманаха в формате YUMA. На CD-диске программа находится в папке 06_Vsion_GLONASS_GPS. Файл альманаха мож- но получить на сайте Национального авиационного университета или запросить по элек- тронной почте cnsatm@nau. edu.ua. Цель лабораторной работы изучение орбитального движения навигационных спутников GPS и ГЛОНАСС с последующей визуализацией спутников на любой момент времени. 5.3.2 Лабораторная работа 5. 3 «Модель движения и визуализация спутников GPS и ГЛОНАСС» Рекомендуется следующий порядок выполнения лабораторной работы. 1. Создайте папку Vsion_GLONASS_GPS_My. 2. Скопируйте в папку Vsion_GLONASS_GPS_My из Vsion_GLONASS_GPS функции и файлы. 173 3. Запишите в папку Vsion_GLONASS_GPS_My альманах спутников GPS и ГЛОНАСС в формате YUMA и выполните задания 1, 2, 3 4. Задание 1. Разберите тексты функций и файлов в папке Vsion_GLONASS_GPS_My. Откройте файл Vision_GLONASS_GPS.m и внесите в него входные данные: имя аль- манаха и время, на которое выполняется моделирование (в комментариях к файлу это указано). Выполните файл и результаты визуализации с полярной диаграммы внесите в отчет. 5. Задание 2. Измените время наблюдения на 6 часов, выполните файл Vision_GLONASS_GPS.m. Занесите результаты наблюдения в отчет и объясните при- чины изменения конфигурации спутников. 6. Задание 3. В файле Vision_GLONASS_GPS.m. измените дату наблюдения на два дня, сохранив время наблюдения и выполните файл. Занесите результаты выполнения про- граммы в отчет и объясните причины, по которым спутники ГЛОНАСС изменили свое местоположения, а спутники GPS почти остались на том же месте. 5.3.3 Контрольные вопросы и задания для самоподготовки 1. Чему равен период обращения спутников GPS?. 2. Чему равен период обращения спутников ГЛОНАСС? 3. Применяя второй закон Кеплера и данные альманаха определите периоды обращения спутников GPS и спутников ГЛОНАСС. 5.3.4 Функции и файлы из папки Vsion_GLONASS_GPS Файл Vsion_GLONASS_GPS.m clear all %Имя m-файла:Vsion_GLONASS_GPS.m %Программа рассчитывает углы видимости навигационных спутников GPS и ГЛОНАСС (азимута и места) %видимых спутников на заданный момент времени %Входные данные: %файл альманаха в формате Yuma,имя файла альманаха присваивается %переменной "Dat",например,Dat = 'имя файла альманаха'; %данные о начале отсчета "d2",d2='месяц/день/год';h=час;min=минута;s=секунда; %координаты позиции приемника -lat(широта в радианах),lon (долгота в радианах, %hr (высота в метрах); %шаг с каким будут рассчитываться параметры (step,секунды); %количество точек (L), в которых будут рассчитываться параметры орбит %L=12*3600/step,L читается так: количество часов (например,12) 174 %число секунд в часе (3600) деленное на шаг (step) %Постоянные: %скорость вращения Земли OMEGAeDOT=7.2921151467e-005; %радиусы земного эллипсоида A_WGS84=6378137.0; B_WGS84=6356752.314; %константы mu=398600500000000; F_CONST = 4.442807633E-10; kt=1;%установка времени на титульной надписи графика, определяется параметрами d2'; h; min; s и j или L; %Задание цветов для графики color6(1:14) = ['k' 'b' 'g' 'r' 'c' 'm' 'r' ':' 'g' ':' 'b' ':' 'k' 'h']; %Входные данные %Dat = 'AlmGLONASS11_12_06.yum'; Dat = 'AlmGG_1.yum'; %Dat = 'YumaGL8_11_06.alm'; d2='12/21/2006'; h=14; min=52.0; s=58.0; %координаты точки, из которой наблюдаются спутники lat = 0.88032730015257;%50 градусов 26минут 20.54 секунд lon = 0.53109641675259;%30 градусов 25 минут 46.4995секунд hr=187.488;% метров X_label=['Широта' ':' num2str(lat) ';' 'долгота' ':' num2str(lon) ';' 'высота' ':' num2str(hr)]; step=300;%шаг отсчета времени в секундах (300=5 минутам);шаг можно изменять step=0; %L=(24*3600) / step;% убрать % перед L для вывода таблицы улов видимости и азимута L=1;% установить % перед L для вывода таблицы улов видимости и азимута %Чтение альманаха [alm,max_kol] = Yuma_GPS_GLONAS_Alm(Dat); %max_kol = 32; nom = 1; i = 0; k = 0; i = 0; % for i = 1 : max_kol while ( i <61) 175 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); nom_ns % - номер навигационного спутника [Rx,Ry,Rz] = ECEFLLH(lon, lat,hr); %Rx=0;Ry=0;,Rz=0;%центр масс Земли %Начало отсчета текущего времени [week,modeweek,d,dweek,weeks]=Tim(d2,h,min,s); for j = 1% 0:L % for j = 1:L % 0:L t( j )=weeks+step*(j); %-step; %t1(j) = t(j)/60; %изменение дискретности текущего времени %d_wn = (week - alm(i).Week); %d_wn = 0; for k = 1 : kol i = nom_ns(k) ; % input "i" !!! % fprintf('i=%i alm(i).Week= %i \n',i, alm(i).Week) d_wn=(modeweek - alm(i).Week);%еслив альманахе не учтено 1024 tk = t(j) + d_wn * 604800 - alm(i).TOA; d_wn = abs(modeweek - alm(i).Week); 176 dd = 302400.0 + d_wn * 604800; if ( ( alm(i).A05 > 0 ) & ( alm(i).Health == 0 ) ) 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)); j2 = 1082.68E-6; re = (A_WGS84 + B_WGS84) / 2.; sin55 = sin(55.0 * pi / 180.0); dn = 1.5 * j2 * re * re / (alm(i).A05^4 ) * (1. - 1.5 * sin55 * sin55); %dn = 0; n=n0 * (1 + dn); 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) 177 break end x = (x2 * y - x * y1) / (y - y1); end % kepler 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; vd = 1. - alm(i).e * cos(Ek); nuk =atan2(sqrt(1-alm(i).e^2)*sin(Ek) / vd,(cos(Ek)-alm(i).e) / vd); Ek = acos((alm(i).e+cos(nuk))/(1+alm(i).e*cos(nuk))); Fk =nuk + alm(i).omega; uk =Fk; ik=alm(i).deltai; % alm(i) 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; % fprintf('i=%i j=%i \n ',i, j); Xk(j,i) = xkk*cos(OMEGAk)-ykk*cos(ik)*sin(OMEGAk); Yk(j,i) = xkk*sin(OMEGAk)+ykk*cos(ik)*cos(OMEGAk); Zk(j,i) = ykk*sin(ik); %Дальности до спутников PR(j,i) = sqrt((Xk(j,i) - Rx)^2 + (Yk(j,i) - Ry)^2 + (Zk(j,i) - Rz)^2); %Перевод в географическую систему если требуется %[lons,lats,hrs] = LLHECEF(Xk,Yk,Zk); %(Llon(j),Llat(j),Hhr(j)) = [lons,lats,hrs]; %расчет угла видимости спутника xls = Xk(j,i) - Rx; yls = Yk(j,i) - Ry; zls = Zk(j,i) - 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; 178 xll = xls/range1; yll = yls/range1; zll = zls/range1; if tdot >= 1.00 b = 0.0; elseif tdot <= -1.00 b = pi; else b = acos( tdot); end satang = pi/2.0 - b; TT =satang; TT(j,i) =TT;%угол видимости спутников %расчет угла азимута спутников 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,i) = atan2(xaz,yaz); end if az(j,i) < 0 az(j,i) = az(j,i) + pi*2; end AZ(j,i) =az(j,i) *180/pi;%угла азимута спутников в градусах EL(j,i) = TT(j,i) *180/pi;%угла видимости спутников в градусах % ПЕРЕСЧЕТ ВРЕМЕНИ A(j)=mod(t(j),86400); her(j)=floor(A(j)/3600); m(j)=floor((A(j)-her(j)*3600)/60); sek(j)=A(j)-her(j)*3600-m(j)*60; %Построение полярной системы координат 179 if EL(j,i) < 0 elp = 180; else elp = (EL(j,i) - 90); end; azp = (AZ(j,i) + 90.0); rad = pi / 180; x0 = 0; y0 = 0; xt(j,i) = (elp * cos(azp * rad)); yt(j,i) = -(elp * sin(azp * rad)); % fprintf('i=%i j=%i \n' , i, j); end % i = ns end; % if ( alm(i).A05 > 0 ) end % j = time %ВНИМАНИЕ. Для вывода времени визуализации спутников на график установите kt t_itle=[d2 ' ' num2str(her(kt)) ':' num2str(m(kt)) ':' num2str(sek(kt))]; %X_label=['Широта' ':' num2str(lat) ';' 'долгота' ':' num2str(lon) ';' 'высота' ':' num2str(hr)]; %num2ctr(lat) %num2str(her(kt)) %X_label=['66' ':']; n = 6; max_n = max(nom_ns); n_end = mod(max(nom_ns),n); n_end = mod(kol, n); n2 = fix(kol / n) * n - n +1; %Формирование таблицы вывода времени UTC (Time), GPS (Tgps в секундах), номера спутника (nns), % углов видимости и азимута от времени и номера спутника for i=1:n:kol fprintf(' Time Tgps '); for k=1: n nns = nom_ns(i+k-1); fprintf(' %2i ', nns); end; fprintf(' \n'); for j=1:L fprintf('%2i:%2i:%2i %i ',her(j),m(j),sek(j), t(j)); for k=1: n nns = nom_ns(i+k-1); fprintf('%6.1f *%6.1f ', EL(j,nns), AZ(j,nns) ); 180 end; fprintf(' \n'); end ; % J=1:L if (i) == (n2) n = n_end; end; end% i hold on %Окружности уровней на круговой диаграмме видимости спутников k1 = 10; k2 = 30; k3 = 50; k4 = 70; k5 = 85; k6=90; n=0; for k=1:5:365 n=n+1; m1 = pi / 180; x(n)=cos(k*m1); y(n)=sin(k*m1); end; %График круговой диаграммы plot(k1*x(:),k1*y(:),'k:', k2*x(:),k2*y(:),'k:', k3*x(:),k3*y(:),'k:',k4*x(:),k4*y(:),'k:', k5*x(:),k5*y(:),'r', k6*x(:),k6*y(:),'r:'); text(5, 10,'80','FontSize',12,'FontName','TimesNewRoman'); text(18, 23,'60','FontSize',12,'FontName','TimesNewRoman'); text(32, 37,'40','FontSize',12,'FontName','TimesNewRoman'); text(45, 50,'20','FontSize',12,'FontName','TimesNewRoman'); text(55, 60,'5','FontSize',12,'FontName','TimesNewRoman'); text(62, 67,'0','FontSize',12,'FontName','TimesNewRoman'); grid on; %Построение изображений видимых спутников на круговой диаграмме i=1; for k=1:kol i = nom_ns(k) ; if (i < 31) plot(xt(kt,i),yt(kt,i), 'Marker' ,'o','MarkerSize',20,'MarkerFaceColor','g' ) ; else plot(xt(kt,i),yt(kt,i), 'Marker' ,'o','MarkerSize',20,'MarkerFaceColor','r' ) ; end; 181 title(t_itle); xlabel(X_label,'FontSize',12,'FontName','TimesNewRoman') set(get(gcf,'CurrentAxes'),'FontSize',14,'FontName','TimesNewRoman') hold on if i<32 str1 = num2str( i, 2); else str1 = num2str( (i - 37), 2); end text(xt(kt,i), yt(kt,i),str1,'FontSize',14,'FontName','TimesNewRoman','HorizontalAlignment','center' ); axis( [-100 100 -100 100]); %axis( [-90 90 -90 90]); end clear Функция LLHECEF function [lons,lats,hrs] = LLHECEF(Xk,Yk,Zk) %Имя функции: LLHECEF %Функция выполняет преобразование координат. %Входные данные:Xk (Rx), Yk (Ry), Zk (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 Функция ECEFLLH function [Rx,Ry,Rz] = ECEFLLH(lon, lat,hr) %Имя функции: ECEFLLH %Функция выполняет преобразование координат 182 %Входные данные:lon-долгота,lat-широта,h-высота;a,b-большая %и малая полуоси эллипсоида %Выходные данные:Rx,Ry,Rz- координаты в ECEF %Для WGS-84 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); Функция Tim function [week,modeweek,d,dweek,weeks]=Tim(d2,h,min,s) %Имя функции:Tim %Функция по данным d2='месяц/день/год';h=час;min=минута;s=секунда рассчитывает % номер недели ([week), при необходимости номер модифицированной недели (modeweek), число дней, %прошедших с 6 января 1980 года (d), номер дня недели (dweek), время GPS в секундах от начала не- дели. %Неделя начинается в ночь с субботы на воскресенье. Воскресенье нулевой день. %d2='10/35/2003';h=23.0;min=59.0;s=59.0; %week = ceil(DAYSDIF('1/6/1980',d2,3)/7); week = floor(DAYSDIF('1/6/1980',d2,3)/7);% текущая неделя modeweek=week; %modeweek=week - 1024; d = DAYSDIF('1/6/1980',d2,3); dweek=fix(d-week*7); weeks=(dweek)*24*60*60+h*60*60+min*60+s; Результат выполнения программы изображен на рис. 5.2. 183 -100 -50 0 50 100 -100 -80 -60 -40 -20 0 20 40 60 80 100 80 60 40 20 5 0 12/21/2006 14:52:58 Широта:0.88033;долгота:0.5311;высота:187.488 3 7 14 16 18 19 21 22 26 29 1 7 8 22 24 Рис. 5.2 . Видимость спутников GPS и ГЛОНАСС На рис. 5.3 показаны экспериментальные данные видимых спутников, полученные приемником НАВИОР 14. Рис. 5.3 . Видимость спутников GPS и ГЛОНАСС (данные измерений) 184 Приведенная программа позволяет также рассчитать углы азимута и места в течение интервала времени, задаваемого пользователем и прогнозировать сеансы навигационных определений при совместном применении систем GPS и ГЛОНАСС. 5.4 Модель движения спутников GPS, ГЛОНАСС, GALILEO 5.4.1 Краткие сведения из теории Рассмотрим кратко общие сведения о системе GALILEO, которая в настоящее время находится в стадии создания. Европейская спутниковая система GALILEO будет содер- жать 27 рабочих и три резервных навигационных спутника, расположенных равномерно на трех орбитах . 5.4 ( рис ) Рис. 5.4. Орбитальное положение спутников GALILEO (модель Долготы восходящих узлов орбит отстоят друг от друга на 120 °. В плоскости орбит спутники расположены через 40 °. Наклонение орбит к плоскости экватора составляет 56°. Средние значения больших полуосей орбит и времени обращения составляют 30049415.54 м и 12 часов 24 мин. соответственно. Система GALILEO будет предоставлять следующие виды обслуживания: открытое обслуживание, OS (Open Service)- сигналы доступные всем пользователям; служба поиска и спасения, SAR (Search and Rescue)- сигналы, доступные всем поль- зователям; служба безопасности движения, SLS (Safety of Life Service) – сигналы доступны авиационным и морским потребителям на договорной основе; коммерческая служба, CS (Commercial Service) – сигналы доступны пользователям на платной основе; ) 185 государственное регулируемое обслуживание, PRC (Public Regulated Service)- го- сударственным службам предоставляются помехозащищенные и зашифрованные навига- ционные сигналы. Данные об орбитальном движении спутников GALILEO будут отличаться от дан- ных GPS только количественными значениями. Поскольку эта система еще не развернута, то в рассматриваемой ниже модели в интегральном альманахе GPS, ГЛОНАСС, GALILEO данные по GALILEO получены искусственно на основании сведений, имеющихся в пе- риодических публикациях. |