实用文档之卫星导航定位算法与程序设计——实验报告

1、实用文档之2013级测绘工程专业卫星导航定位算法与程序设计实验报告实验名称：卫星导航基本程序设计班 级：学 号：姓 名：实验时间：2016年6月28日2016年6月30中国矿业大学目录实验一时空基准转换2一、实验目的2二、实验内容2三、实验过程2四、实验感想6实验二 RINEX文件读写7一、实验目的7、实验内容7三、实验过程8实验三卫星轨道计算13一、实验目的13二、实验内容13三、实验过程13四、实验感想16实验一时空基准转换一、实验目的1、加深对时空系统及其之间转换关系的理解2、掌握常用时空基准之间的转换模型与软件实现3、每人独立完成实验规定的内容二、实验内容本实验内容包括：内容一：编程实

2、现GPS起点1980年1月6日0时对应的儒略日 内容二：编程实现2011年11月27日对应的GPS周数与一周内 的秒数内容三：在WGS84椭球的条件下，编程实现当中央子午线为117 度时，计算高斯坐标 x = 3548910.811290287, y = 179854.6172135982 对应的经纬度坐标？内容四：WGS84 椭球下，表面 x=-2408000; y=4698000;z= 3566000 处的地平坐标系坐标为：e=704.8615;n=l 14.8683;u=751.9771的点对应 的直角坐标为多少？三、实验过程1针对第一、二部分内容：1.1解决思路:先建立” TimeSt

3、ruct.h”的头文件,将格里高利历、GPS时间结构、儒略日时间结构共结构体的方式放在里面；在建立 “TimeTr”的头文件，建立类“CTinwTr”，创建变量“GPSTime”、 “Time”、”JulDay”，并且申明函数“TIME2JUL”、“TIME2GTIME”等， 用这些函数分别实现所需要的转换。1.2具体的实现函数:“TIME2JUL” 函数：double CTimeTr:TIME2JUL 0/TINE Time, JULIAXDAY &JulDaydouble m, y;double D;/h 二Time. byHour+Time. byMinute/60. O+Time. d

4、Second/3600. 00;if(Time. byMonth=2)y=Time wYearl;m=Time. byMonth+12;elsey=Time wYear;m=Time byMonth:D=floor(365. 25*(y+4716)+floor(30. 6001*(m+1)+Time. byDay+Time. byHo ur/24. 0-1537. 5;JulDay. IDay = int(D);JulDaylSecond = Dint(julDayIDay);return 0;“TIME2GTIME”：void CTimeTr: :TIME2GTIME()double JD;

5、long m, y;int WN;double Wsecend;/UT=Time. byHour+Time. byMinute/60. 0+Time. dSecond/3600. 00;if(Time. byMonth=00001 & abs (ll-10)=0.0001BO=B1;aa0=aal;10=11;F_:B= (c*Beta2+ (c*Beta4+ (c*Be,ta6+c*Be,ta8*cos (BO) A2) *c os (BO)A2)*cos (BO)A2)*sin(BO)*cos(BO);l=a2*10A2+a4*10A4+a6*10A6;F_yBl=a3*10A3+a5*1

6、0A5;Bl=(x-F_xB-F_xBl)/(c*Beta0);aal=(a*cos(Bl)/sqrt(l-eA2*(sin(Bl)A2); 11=(y-F_yBl)/aal;endL=rad2deg(11)+L0;B=rad2deg(Bl);2.3实验结果B L Mffl Lxl double1234567132.04922B L士 lxl double12345671118.90422四、实验感想本次试验是花时间较多的一次实验，关于时间转换的部分全部都 是自己动手将matlab代码写成“C+”的类，进行实现的。其中遇到的 较大的困难是儒略日向UTC转换的部分，这部分的函数步骤较多， 关键是

7、在一开始的时间结构里面，各时间各部分的数据类型大多定义 的是“int”型的，但是在进行计算的时候有较多的小数，需要用到浮点 型的函数，这部分用了较多的时间。在做这个实验的时候，第一天花了时间主要是转换代码，使程序 没有错误，能够正常的运行出来，出现黑框框，但是还只有个别功能 能够用，能够运行出正确的结果：第二天时间主要是花在修改函数上 头，能够使所写的功能都能运行出正确的结果。通过做时间转换的实 验，使自己产生了第一次亲自编写“C+”代码的经历，而且所有错误 的解决全部都是自己解决，收获不少。实验二RINEX文件读写一、实验目的1、深入了解RINEX文件格式2、进一步提高MATLAB程序设计能

8、力3、掌握N文件、O文件、SP3文件的基本读写技巧二、实验内容本实验内容包括：1、任选IGS站，下载N文件、O文件与SP3文件；2、编程实现N文件读入，并采用中文标注出主要参数的名称及 作用；4、编程实现O文件读入，并釆用中文标注出主要参数的名称及作用:5、编程实现SP3文件读入，并釆用中文标注出主要参数的名称 及作用；三、实验过程1、针对第一部分内容：编程实现N文件读入，并采用中文标注出主要参数的名称及作用1. 1、解决思路：按照“GPSeasy开源代码提供的函数，按照实验要求读取了 N文件 的内容，先用“rinexe”函数，将N文件读取成“eph.dat”文件，然后再 用Mget_eph函

9、数将“eph. dat”文件读取成“Eph”矩阵，此矩阵中包含了 N 文件中的数据，在最后用“fprintf”函数将所需要的数据输出成TXT 文件即可。1.2、主要函数代码：“get_eph” 函数：function eph = get_eph(ephemeridesfile)fide = fopen(ephemeridesfile);eph, count = fread(fide, Inf, 1 double 1);noeph = count/22;eph = reshape (ephz 22f noeph);“rinexe”函数(部分)：function rinexe(ephemerisf

10、i1e, outputfile)fide = fopen(ephemerisfile);headlines = 0;while 1head_lines = head_lines+l;line = fgetl(fide);answer = finds trains END OF HEADER);if -isempty(answer), break; end;end;headlines主函数中输出结果得函数部分:afO=data (19)川卫星中差MO=data(3);roota=data(4);deltan=data(5);ecc=data(6);omega=data(7);cuc=data (

11、8);cus=data(9);crc=data(10);crs=data(11);i0=data(12);idot=data(13);toe=data(18);afl=data (20);务对所要输出的参数赋值fprintf (fid, 1 n卫星编号：%dn卫星钟差：%dn平近点角距：%dn轨 道长半轴的平方根：%dn平均运动修正量：%dn轨道偏心率：%dn近地点角 距：%dn 纬度幅角的余弦调和项改正的振幅 1,prn,af0,M0z roota,deltan，ecc,omega,cue);fprintf (fid,，纬度幅角的正弦调和项改正的振幅：%dn轨道半径的余弦 调和项改正的振幅：

12、%dn轨道半径的正弦调和项改正的振幅：%dn轨道倾 角：%dn 轨道倾角变化率：%dn 星历参考时 刻:%dn ,cus,crc,ers,iO,idotz toe)fclose(fid);1.3.输出结果2、针对第二部分内容：编程实现O文件读入，并釆用中文标注出主要参数的名称及作 用；2.1、实现思路：通过matlab的函数“fopen”读取O文件，得到O 文件的指针，通过“anheader”函数将文件中的接收机大致位置” approx_XYZlv ,天线的偏移值ant_deltal,观测值类型 obstypesl等读入成为matlab的矩阵，然后通过循环，利用 “grabdata”函数将所需

13、要的历元的观测文件依次输出来，最后通过函数，将所需要的数据依次打印出来。2.2、主要函数：anheader” 函数：funotion Obs_types, ant_delta# ifound_types,approx_XYZ = anheader(file)fid = fopen (file, 1rt1);eof = 0;ifound_types = 0;Obs_types =;ant_delta =; approxXYZ =;while 1% Gobbling the headerline = fgetl(fid);answer = findstr(line, 1 END OF HEADER

14、1);if -isempty (ansvzer) , break; end;if (line = -1) , eof = 1; break; end;answer = findstr(line,ANTENNA: DELTA H/E/N1);if isempty(answer)for k = 1:3delta, line = strtok(line);del = str2num(delta);ant_delta = ant_delta del;end;endanswer = findstr(line,APPROX POSITION XYZ1);if isempty(answer)for k =

15、1:3app_XYZz line = strtok (line);del = str2num(app_XYZ); approx_XYZ = approx_XYZ del;end;endanswer = findstr(line,# / TYPES OF OBSERV1);if isempty(answer)NObs r line = strtok(line);NoObs = str2num(NObs);for k = 1:NoObsot, line = strtok(line);Obs_types = Obs_types ot;end;ifound_types = 1;end;end;%fcl

16、ose(fid);grabdata” 函数：function Obs = grabdaNoSv, NoObs)%GRABDATA Positioned in a RINEX file at a selected epoch%reads observations of NoSv satellitmsglobal linObs = zeros(NoSvz NoObs);if NoObs = 5% This will typical be Turbo SII datafor u = 1:NoSvlin = fgetl(fid);for k = 1:NoObsObs(uzk) = str2num(li

17、n(2+16*(k-1):16*k-2);endendelse% This will typical be Z12 dataObs = Obs (: , 1 2 3 4 5 ) ; % We cancel rhe last mo columns6 and 7NoObs = 5;for u = 1:NoSvlin = fgetl(fid);lin_doppler = fgetl(fid);for k = 1:NoObs %-lif isempty(str2num(lin(l+16*(k-1):16*k-2) = 1,Obs(u,k) = nan;else%Obs(u,k) = str2num(l

18、in(1 + 16*(k-1) :16*k-2); end% Obs(u,NoObs) = str2num(lin(65:78);endendend2.3实验结果四、实验感想这部分实验是我在之前做的，之前自己有看过“gps_easy”有关的 代码，看过相关的“N文件” “O文件”读写函数，并且学会了如何 调用这些函数，对里面的输出量有了一点的了解，所以我自己的主要 工作就是运用了“rintF函数，将读取到matlab中的矩阵写入TXT文 档中，这部分工作量不是很大，但较有意义。实验三卫星轨道计算一、实验目的1、进一步熟悉N文件的读入2、掌握开普勒参数讣算卫星轨道的过程3、编程实现采用广播星历计

19、算卫星轨道4、掌握MATLAB函数调用步骤二、实验内容本实验内容包括：1、调试时间转换函数，熟悉内容，备主函数调用2、调试广播星历导航文件的读入程序，备主函数调用3、根据卫星位置计算公式编写主函数，同时调用时间转换、星 历读取等的子函数来共同完成卫星位置的计算，最后输出讣算结果4、理清程疗:各模块的功能结构三、实验过程1、实验思路：在老师提供的“SPP”文件中，直接利用卫星位置计算函数，进行 卫星位置的计算，将利川Gps.cpp文件中的GetGpsPosition函数，利 用其中的迭代求解卫星位置部分，用“cout”直接将卫星迭代后的位置 直接输出，因为星历文件中有较多的星历，所以利用循环语句

20、，将求 解出来的卫星位置依次输出出来。2、主要函数bool CGps：:GetGpsPosition()GPSTIME ts;GPSTIME tr;GPSTIME tsO;GPSTIME oGTime;GPSTIME nGTOC;vectorSendSignPosition;GpsPos GpsPTemp;GpsSendPosition SdSignPoTemp;int nTheFitPoint=0;Position pTemp;GpsPosition clear 0;if (oData. size () =0)strErr = (PRN=%d没有对应星历)；return false:if (

21、nData. size 0=0)return false:coutfixed;for (int i=0;ioDatasize();i+)double *VX=new double ；/for (int j=0;jle-007);/earthrot(trlSecond-tsISecond, pTemp XX, pTemp YY, pTemp ZZ);TimeToGpsTime(nDataZnTheFitPoint. TOC, nGTOC);/SdSignPoTemp dt=nDatanTheFitPointZ dClkBias+nData InTheFitPoint. dClkDrift*(ts

22、. lSecondnGTOC. ISecond) + nData InTheFitPoint dClkDriftRate*pow(ts. ISecond-nGTOC. ISecond),2);/SdSignPoTemp. PRN=oDatai. oObserveDataj. PRN;SdSignPoTemp XX=pTemp XX;SdSignPoTemp. YY=pTemp. YY;SdSignPoTemp ZZ=pTemp ZZ;SdSignPoTemp SendTime lWeek=oGTinie lWeek;SdSignPoTemp SendTimeISecond二tsISecond;

23、SdSignPoTemp CA=oDatai oObserveDataj dCl;if(i10)couti+l z z，SdSignPoTempSdSignPoTemp. XX*SdSignPoTemp. YY2 16135223.755143 185584.575777 143? 371B3V .71134? 5859516.861784 16349120.077991 19947316.778837 1439 371039.911790 -2R784709 42642 9125024.1450R8 14R471231439 371039.91 097R16155592.G27192 209

24、10315.932008 -2089955.58830S 1439 3710549i.lG1.53 -7247607.767861 2449?492.938252 -3356507.733811 1439 371054.?119487 -1.5276249.624856 -5068117.335053 21478068.438348 1439 371054.91048516 -7226927.388777 161fi255H.222SB7 1987S897.371429 1439 371354.91134623 5837736.943839 16380673.518278 19927561.864706 1439 371054.91179031 -20810783.607651 9121236950丄27 14010542.604511 1439 371054.91.D9691 -1615555U.6b869 2MVB516M.8b6kr/6 -20S8428.35?7S9 1439 371W69.911615 3 -9255870.265311 24502184.913283 -3310321.9