卫星导航定位算法与程序设计实验报告

2013 级测绘工程专业卫星导航定位算法与程序设计实验报告实验名称：卫星导航基本程序设计班级：学号：姓名：实验时间： 2016年6月28日~2016年6月30中国矿业大学目录实验一时空基准转换 (2)一、实验目的 (2)二、实验内容 (2)三、实验过程 (2)四、实验感想 (6)实验二RINEX文件读写 (6)一、实验目的 (6)二、实验内容 (7)三、实验过程 (7)实验三卫星轨道计算 (12)一、实验目的 (12)二、实验内容 (12)三、实验过程 (12)四、实验感想 (15)实验一时空基准转换一、实验目的1、加深对时空系统及其之间转换关系的理解2、掌握常用时空基准之间的转换模型与软件实现3、每人独立完成实验规定的内容二、实验内容本实验内容包括：内容一：编程实现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解决思路：先建立” TimeStruct.h”的头文件，将格里高利历、GPS时间结构、儒略日时间结构共结构体的方式放在里面；在建立“TimeTr”的头文件，建立类“CTimeT r”，创建变量“GPS Time”、“Time”、”JulDay”，并且申明函数“TIME2JUL”、“TIME2GTIME”等，用这些函数分别实现所需要的转换。

1.2具体的实现函数：“TIME2JUL”函数：double CTimeTr::TIME2JUL()//TIME Time,JULIANDAY &JulDay{double m,y;double D;//h =Time.byHour+Time.byMinute/60.0+Time.dSecond/3600.00;if(Time.byMonth<=2){y=Time.wYear-1;m=Time.byMonth+12;}else{y=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 WN;double Wsecend;//UT=Time.byHour+Time.byMinute/60.0+Time.dSecond/3600.00;if(Time.byMonth<=2){y=Time.wYear-1;m=Time.byMonth+12;}else{y=Time.wYear;m=Time.byMonth;}JD=int(365.25*y)+int(30.6001*(m+1))+Time.byDay+Time.byHour/24.0+1720981.5;WN = floor((JD-2444244.5)/7.0);GpsTime.lWeek=WN;Wsecend=(JD-2444244.5-7*WN)*604800;GpsTime.lSecond=Wsecend;}1.3实验结果：2 针对第三部分内容：2.1解决思路：运用实验指导书中提供的matlab高斯反算的代码，进行解算；将高斯反算的公式直接输成matlab代码，绕后在函数“function [B,L] = gauss_fansuan（x,y,L0）”中，将坐标x = 3548910.811290287，y = 179854.6172135982，L0 = 117，带入函数的坐边，即可得到所需要的经纬度。

2.2主要函数的代码：function [B,L]=gauss_fansuan(x,y,L0)a=6378137;f=1/298.257223563;b=a-a*f;c=a^2/b;e=sqrt(a^2-b^2)/a;e1=sqrt(a^2-b^2)/b;Beta0=1-(3/4)*e1^2+(45/64)*e1^4-(175/256)*e1^6+(11025/16384)*e1^8;Beta2=Beta0-1;Beta4=(15/32)*e1^4-(175/384)*e1^6+(3675/8192)*e1^8;Beta6=-(35/96)*e1^6+(735/2048)*e1^8;Beta8=(315/1024)*e1^8;B0=x/(c*Beta0);aa0=(a*cos(B0))/sqrt(1-e^2*sin(B0)^2);l0=y/aa0;N=a*sqrt(1-e^2*sin(B0)^2);t=tan(B0);in=e1*cos(B0);a1=N*cos(B0);a2=(1/2)*N*sin(B0)*cos(B0);a3=(1/6)*N*cos(B0)^3*(1-t^2+in^2);a4=(1/24)*N*sin(B0)*cos(B0)^3*(5-t^2+9*in^2+4*in^4);a5=(1/120)*N*cos(B0)^5*(5-18*t^2+t^4+14*in^2-58*in^2*t^2);a6=(1/720)*N*sin(B0)*cos(B0)^5*(61-58*t^2+t^4);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*l0^2+a4*l0^4+a6*l0^6;F_yBl=a3*l0^3+a5*l0^5;B1=(x-F_xB-F_xBl)/(c*Beta0);aa1=(a*cos(B1))/sqrt(1-e^2*(sin(B1))^2);l1=(y-F_yBl)/aa1;while abs(B1-B0)>=0.0001 && abs(l1-l0)>=0.0001B0=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*l0^2+a4*l0^4+a6*l0^6;F_yBl=a3*l0^3+a5*l0^5;B1=(x-F_xB-F_xBl)/(c*Beta0);aa1=(a*cos(B1))/sqrt(1-e^2*(sin(B1)^2));l1=(y-F_yBl)/aa1;endL=rad2deg(l1)+L0;B=rad2deg(B1);2.3实验结果四、实验感想本次试验是花时间较多的一次实验，关于时间转换的部分全部都是自己动手将matlab代码写成“C++”的类，进行实现的。

其中遇到的较大的困难是儒略日向UTC转换的部分，这部分的函数步骤较多，关键是在一开始的时间结构里面，各时间各部分的数据类型大多定义的是“int”型的，但是在进行计算的时候有较多的小数，需要用到浮点型的函数，这部分用了较多的时间。

在做这个实验的时候，第一天花了时间主要是转换代码，使程序没有错误，能够正常的运行出来，出现黑框框，但是还只有个别功能能够用，能够运行出正确的结果；第二天时间主要是花在修改函数上头，能够使所写的功能都能运行出正确的结果。

通过做时间转换的实验，使自己产生了第一次亲自编写“C++”代码的经历，而且所有错误的解决全部都是自己解决，收获不少。

实验二RINEX文件读写一、实验目的1、深入了解RINEX文件格式2、进一步提高MATLAB程序设计能力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”文件，然后再用“get_eph”函数将“eph.dat”文件读取成“Eph”矩阵，此矩阵中包含了N文件中的数据，在最后用“fprintf”函数将所需要的数据输出成”.TXT”文件即可。

1.2、主要函数代码:“get_eph”函数：function 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 1head_lines = head_lines+1;line = fgetl(fide);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);%对所要输出的参数赋值fprintf(fid,'\n卫星编号：%d\n卫星钟差：%d\n平近点角距：%d\n轨道长半轴的平方根：%d\n平均运动修正量：%d\n轨道偏心率：%d\n近地点角距：%d\n纬度幅角的余弦调和项改正的振幅',prn,af0,M0,roota,deltan,ecc,omega,cuc);fprintf(fid,'纬度幅角的正弦调和项改正的振幅：%d\n轨道半径的余弦调和项改正的振幅：%d\n轨道半径的正弦调和项改正的振幅：%d\n轨道倾角：%d\n轨道倾角变化率：%d\n 星历参考时刻：%d\n',cus,crc,crs,i0,idot,toe)fclose(fid);1.3、输出结果2、针对第二部分内容：编程实现O文件读入，并采用中文标注出主要参数的名称及作用；2.1、实现思路：通过matlab的函数“fopen”读取O文件，得到O文件的指针，通过“anheader”函数将文件中的接收机大致位置”approx_XYZ1”，天线的偏移值”ant_delta1”，观测值类型“Obs_types1”等读入成为matlab的矩阵，然后通过循环，利用“grabdata”函数将所需要的历元的观测文件依次输出来，最后通过“fprintf”函数，将所需要的数据依次打印出来。