《卫星导航定位算法与程序设计实验报告.doc》由会员分享,可在线阅读,更多相关《卫星导航定位算法与程序设计实验报告.doc(20页珍藏版)》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、.-2013 级测绘工程专业卫星导航定位算法与程序设计实验报告实验名称: 卫星导航基本程序设计 班 级: 学 号: 姓 名: 实验时间: 2016年6月28日2016年6月30 中 国 矿 业 大 学目录实验一 时空基准转换2一、实验目的2二、实验内容2三、实验过程2四、实验感想6实验二 RINEX文件读写7一、实验目的7二、实验内容7三、实验过程7实验三 卫星轨道计算12一、实验目的12二、实验内容12三、实验过程12四、实验感想15实验一 时空基准转换一、实验目的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=114.8683;u=751.9771的点对应的直角坐标为多少?三、实验过程1.针对第一、二部分内容:1.1解决思路:先建立”
3、 TimeStruct.h”的头文件,将格里高利历、GPS时间结构、儒略日时间结构共结构体的方式放在里面;在建立“TimeTr”的头文件,建立类“CTimeTr”,创建变量“GPSTime”、“Time”、”JulDay”,并且申明函数“TIME2JUL”、“TIME2GTIME”等,用这些函数分别实现所需要的转换。1.2具体的实现函数:“TIME2JUL”函数:double CTimeTr:TIME2JUL()/TIME Time,JULIANDAY &JulDaydouble m,y;double D;/h =Time.byHour+Time.byMinute/60.0+Time.dSec
4、ond/3600.00;if(Time.byMonth=2)y=Time.wYear-1;m=Time.byMonth+12;elsey=Time.wYear;m=Time.byMonth;D=floor(365.25*(y+4716)+floor(30.6001*(m+1)+Time.byDay+Time.byHour/24.0-1537.5;JulDay.lDay = int(D);JulDay.lSecond = D-int(JulDay.lDay);return 0;“TIME2GTIME”:void CTimeTr:TIME2GTIME()double JD;long m,y;int
5、 WN;double Wsecend;/UT=Time.byHour+Time.byMinute/60.0+Time.dSecond/3600.00;if(Time.byMonth=0.0001 & abs(l1-l0)=0.0001 B0=B1; aa0=aa1; l0=l1; F_xB=(c*Beta2+(c*Beta4+(c*Beta6+c*Beta8*cos(B0)2)*cos(B0)2)*cos(B0)2)*sin(B0)*cos(B0); F_xBl=a2*l02+a4*l04+a6*l06; F_yBl=a3*l03+a5*l05; B1=(x-F_xB-F_xBl)/(c*Be
6、ta0); aa1=(a*cos(B1)/sqrt(1-e2*(sin(B1)2); l1=(y-F_yBl)/aa1;endL=rad2deg(l1)+L0;B=rad2deg(B1);2.3实验结果四、实验感想本次试验是花时间较多的一次实验,关于时间转换的部分全部都是自己动手将matlab代码写成“C+”的类,进行实现的。其中遇到的较大的困难是儒略日向UTC转换的部分,这部分的函数步骤较多,关键是在一开始的时间结构里面,各时间各部分的数据类型大多定义的是“int”型的,但是在进行计算的时候有较多的小数,需要用到浮点型的函数,这部分用了较多的时间。在做这个实验的时候,第一天花了时间主要是转换
7、代码,使程序没有错误,能够正常的运行出来,出现黑框框,但是还只有个别功能能够用,能够运行出正确的结果;第二天时间主要是花在修改函数上头,能够使所写的功能都能运行出正确的结果。通过做时间转换的实验,使自己产生了第一次亲自编写“C+”代码的经历,而且所有错误的解决全部都是自己解决,收获不少。实验二 RINEX文件读写一、实验目的1、深入了解RINEX文件格式2、进一步提高MATLAB程序设计能力3、掌握N文件、O文件、SP3文件的基本读写技巧二、实验内容本实验内容包括:1、任选IGS站,下载N文件、O文件与SP3文件;2、编程实现N文件读入,并采用中文标注出主要参数的名称及作用;4、 编程实现O文
8、件读入,并采用中文标注出主要参数的名称及作用;5、编程实现SP3文件读入,并采用中文标注出主要参数的名称及作用;三、实验过程1、针对第一部分内容:编程实现N文件读入,并采用中文标注出主要参数的名称及作用1.1、解决思路:按照“GPSeasy”开源代码提供的函数,按照实验要求读取了N文件的内容,先用“rinexe”函数,将N文件读取成“eph.dat”文件,然后再用“get_eph”函数将“eph.dat”文件读取成“Eph”矩阵,此矩阵中包含了N文件中的数据,在最后用“fprintf”函数将所需要的数据输出成”.TXT”文件即可。1.2、主要函数代码:“get_eph”函数:function
9、eph = get_eph(ephemeridesfile)fide = fopen(ephemeridesfile);eph, count = fread(fide, Inf, double);noeph = count/22;eph = reshape(eph, 22, noeph);“rinexe”函数(部分):function rinexe(ephemerisfile, outputfile)fide = fopen(ephemerisfile);head_lines = 0;while 1 head_lines = head_lines+1; line = fgetl(fide);
10、answer = findstr(line,END OF HEADER); if isempty(answer), break; end;end;head_lines主函数中输出结果得函数部分: af0=data(19);%卫星中差 M0=data(3); roota=data(4); deltan=data(5); ecc=data(6); omega=data(7); cuc=data(8); cus=data(9); crc=data(10); crs=data(11); i0=data(12); idot=data(13); toe=data(18); af1=data(20);%对所
11、要输出的参数赋值 fprintf(fid,n卫星编号:%dn卫星钟差:%dn平近点角距:%dn轨道长半轴的平方根:%dn平均运动修正量:%dn轨道偏心率:%dn近地点角距:%dn纬度幅角的余弦调和项改正的振幅,prn,af0,M0,roota,deltan,ecc,omega,cuc);fprintf(fid,纬度幅角的正弦调和项改正的振幅:%dn轨道半径的余弦调和项改正的振幅:%dn轨道半径的正弦调和项改正的振幅:%dn轨道倾角:%dn轨道倾角变化率:%dn星历参考时刻:%dn,cus,crc,crs,i0,idot,toe)fclose(fid);1.3、输出结果2、针对第二部分内容:编程
12、实现O文件读入,并采用中文标注出主要参数的名称及作用;2.1、实现思路:通过matlab的函数“fopen”读取O文件,得到O文件的指针,通过“anheader”函数将文件中的接收机大致位置” approx_XYZ1”,天线的偏移值” ant_delta1”,观测值类型“Obs_types1”等读入成为matlab的矩阵,然后通过循环,利用“grabdata”函数将所需要的历元的观测文件依次输出来,最后通过“fprintf”函数,将所需要的数据依次打印出来。2.2、主要函数:“anheader”函数:function Obs_types, ant_delta,ifound_types,appr
13、ox_XYZ = anheader(file)fid = fopen(file,rt);eof = 0;ifound_types = 0;Obs_types = ;ant_delta = ;approx_XYZ = ;while 1 % Gobbling the header line = fgetl(fid); answer = findstr(line,END OF HEADER); if isempty(answer), break; end; if (line = -1), eof = 1; break; end; answer = findstr(line,ANTENNA: DELT
14、A H/E/N); if isempty(answer) for k = 1:3 delta, line = strtok(line); del = str2num(delta); ant_delta = ant_delta del; end; end answer = findstr(line,APPROX POSITION XYZ); if isempty(answer) for k = 1:3 app_XYZ, line = strtok(line); del = str2num(app_XYZ); approx_XYZ = approx_XYZ del; end; end answer
15、 = findstr(line,# / TYPES OF OBSERV); if isempty(answer) NObs, line = strtok(line); NoObs = str2num(NObs); for k = 1:NoObs ot, line = strtok(line); Obs_types = Obs_types ot; end; ifound_types = 1; end;end;%fclose(fid);“grabdata”函数:function Obs = grabdata(fid, NoSv, NoObs)%GRABDATA Positioned in a RI
16、NEX file at a selected epoch% reads observations of NoSv satellitesglobal linObs = zeros(NoSv, NoObs);if NoObs = 5 % This will typical be Turbo SII data for u = 1:NoSv lin = fgetl(fid); for k = 1:NoObs Obs(u,k) = str2num(lin(2+16*(k-1):16*k-2); end endelse % This will typical be Z12 data Obs = Obs(:
17、,1 2 3 4 5); % We cancel the last two columns 6 and 7 NoObs = 5; for u = 1:NoSv lin = fgetl(fid); lin_doppler = fgetl(fid); for k = 1:NoObs %-1 if isempty(str2num(lin(1+16*(k-1):16*k-2) = 1, Obs(u,k) = nan; else % Obs(u,k) = str2num(lin(1+16*(k-1):16*k-2); end % Obs(u,NoObs) = str2num(lin(65:78); en
18、dendend2.3实验结果四、实验感想这部分实验是我在之前做的,之前自己有看过“gps_easy”有关的代码,看过相关的“N文件”“O文件”读写函数,并且学会了如何调用这些函数,对里面的输出量有了一点的了解,所以我自己的主要工作就是运用了“fprintf”函数,将读取到matlab中的矩阵写入TXT文档中,这部分工作量不是很大,但较有意义。实验三 卫星轨道计算一、实验目的1、进一步熟悉N文件的读入2、掌握开普勒参数计算卫星轨道的过程3、编程实现采用广播星历计算卫星轨道4、掌握MATLAB函数调用步骤二、实验内容本实验内容包括:1、调试时间转换函数,熟悉内容,备主函数调用2、调试广播星历导航文
19、件的读入程序,备主函数调用3、根据卫星位置计算公式编写主函数,同时调用时间转换、星历读取等的子函数来共同完成卫星位置的计算,最后输出计算结果4、理清程序各模块的功能结构三、实验过程1、实验思路:在老师提供的“SPP”文件中,直接利用卫星位置计算函数,进行卫星位置的计算,将利用“Gps.cpp”文件中的” GetGpsPosition”函数,利用其中的迭代求解卫星位置部分,用“cout”直接将卫星迭代后的位置直接输出,因为星历文件中有较多的星历,所以利用循环语句,将求解出来的卫星位置依次输出出来。2、主要函数bool CGps:GetGpsPosition()GPSTIME ts;GPSTIME
20、 tr;GPSTIME ts0;GPSTIME oGTime;GPSTIME nGTOC;vectorSendSignPosition;GpsPos GpsPTemp;GpsSendPosition SdSignPoTemp;int nTheFitPoint=0;Position pTemp;GpsPosition.clear();if (oData.size()=0)strErr = (PRN=%d没有对应星历);return false;if (nData.size()=0)return false;coutfixed;for (int i=0;ioData.size();i+)doubl
21、e *VX=new double4;/for (int j=0;j1e-007);/earthrot(tr.lSecond-ts.lSecond,pTemp.XX,pTemp.YY,pTemp.ZZ);TimeToGpsTime(nDatanTheFitPoint.TOC,nGTOC);/SdSignPoTemp.dt=nDatanTheFitPoint.dClkBias+nDatanTheFitPoint.dClkDrift*(ts.lSecond-nGTOC.lSecond)+nDatanTheFitPoint.dClkDriftRate*pow(ts.lSecond-nGTOC.lSec
22、ond),2);/SdSignPoTemp.PRN=oDatai.oObserveDataj.PRN;SdSignPoTemp.XX=pTemp.XX;SdSignPoTemp.YY=pTemp.YY;SdSignPoTemp.ZZ=pTemp.ZZ;SdSignPoTemp.SendTime.lWeek=oGTime.lWeek;SdSignPoTemp.SendTime.lSecond=ts.lSecond;SdSignPoTemp.CA=oDatai.oObserveDataj.dC1;if(i10)couti+1 SdSignPoTemp.PRN SdSignPoTemp.XX SdS
23、ignPoTemp.YY SdSignPoTemp.ZZ SdSignPoTemp.SendTime.lWeek SdSignPoTemp.SendTime.lSecondendl;/SendSignPosition.push_back(SdSignPoTemp);return true;3、实验结果四、实验感想计算卫星位置的这部分“C+”代码比较综合,读代码的难度较大,它综合了前面的“C+”O文件读取,N文件读取,然后后面的计算代码较为复杂,需要对卫星轨道计算的公式了解全面,因此对这部分实验,是在老师和同学的帮助下完成的,直接在轨道计算的函数下加入了输出函数,将计算的轨道数据直接输出出来。通过这部分实验是我加深了对卫星轨道计算公式的了解,以及以后面对复杂的公式应该如何应对。