卫星导航定位算法与程序设计
实验报告
实验一 时空基准转换
一、实验目的
1、加深对时空系统及其之间转换关系的理解
2、掌握常用时空基准之间的转换模型与软件实现
二、实验内容
1、编程实现GPS起点1980年1月6日0时对应的儒略日
2、编程实现2011年11月27日对应的GPS周数与一周内的秒数
3、在WGS84椭球的条件下，编程实现当中央子午线为117度时，计算高斯坐标
x=3548910.811290287，y=179854.6172135982对应的经纬度?
4、在WGS84椭球的条件下, 表面x=-2408000,y=4698000,z=3566000处地平坐标系坐
标为:e=704.8615,n=114.8683,u=751.9771的点对应的直角坐标?
三、实验过程
1、这是测试 站心坐标系 空间直角坐标系 相互转换高斯正算和高斯反算 的主程
序功能模块
（1）调用了enu2xyz 站心坐标向空间坐标的转换函数
（2）调用了xyz2enu空间坐标向站心坐标转换的函数
（3） 调用了高斯反算函数gau_fansuan
2、时间转换主程序功能模块
disp(\'1980年1月6日0时对应的儒略时:\');
jd1 = julday(1980,1,6,0);
disp(jd1);
disp(\'2011年11月27日对应的GPS周和周内秒:\');
jd2 = julday(2011,11,27,0);
[week,sow] = gps_time(jd2);
disp(\'week:\')
disp(week);
disp(\'sow:\');
disp(sow);
四、实验体会
本次实验中我对坐标转换问题有了更多的了解，同时对matlab的一些基本函数更加了解，对时间转换内容的熟悉程度也有所加深。通过本次实验我学会了matlab程序的调试方法，为后继学习打下了基础。
实验二 RINEX文件读写
一、实验目的
1、深入了解RINEX文件格式
2、进一步提高MATLAB程序设计能力
3、掌握N文件、O文件、SP3文件的基本读写技巧
二、实验内容
1、任选IGS站，下载N文件、O文件与SP3文件
2、编程实现N文件的读入，并采用中文标注出主要参数的名称及作用
3、编程实现O文件的读入，并采用中文标注出主要参数的名称及作用
4、编程实现SP3文件的读入，并采用中文标注出主要参数的名称及作用
三、实验过程
1、N文件读取主程序模块功能
（1） 调用rinexe函数将文件SITE247J.01N中的数据按一定的格式读入到文件
eph.dat中
rinexe(\'SITE247J.01N\',\'eph.dat\');
（2） 调用get_eph函数将文件eph.dat中的数据读入到变量Eph中
Eph = get_eph(\'eph.dat\');
（3） 输出得到的Eph变量的内容到命令窗口
disp(Eph);
（4） rinexe函数及get_eph函数具体注释见源代码
2、SP3文件、O文件读取主程序模块功能
读取O文件主要过程：
（1） 打开文件，获取头文件的信息。利用anheader获取观测值类型、天线高及偏心及找到相对应观测类型的标识码和判断文件是否结束的标识码。
（2） 在函数体中利用findstr()及循环语句逐行查找字符串“END OF HEADER”，找到头文件结束的位置。
（3然后逐行在头文件中寻找字符串“ANTENNA: DELTA H/E/N”、”# / TYPES OF OBSERV”等，得到天线高及偏心高矩阵、观测值类型、观测值个数和观测类型的标识码。
（4） 获取观测历元的信息。利用fepoch_0获得某一个历元的观测时间和卫星编号矩阵及文件结束标识符。然后可以从卫星的编号矩阵sats中求得卫星数。
（5） 获取某一观测类型（如程序中的P2）的信息。利用grabdata.m函数,根据卫星数NoSv和观测类型数NoObs返回某历元的5*5\"obs\"观测数据矩阵,再利用循环获得所有历元的观测数据矩阵。
（6） 利用fobs_typ.m获取观测类型为P2的观测值矩阵中对应列的信息。
涉及的主要函数及功能：
anheader.m：分析RINEX文件的头文件，输出观测类型和天线设置类型，及找到相应观测类型的标识码。
fepoch_0.m:在打开的RINEX文件中利用识别标志fid寻找下一历元。从历元序列中生成历元时间、卫星编号、文件结束符。只有O文件才会被处理。
grabdata.m：在选定的历元读取编号卫星对应的观测值。
fobs_typ.m：返回观测值矩阵中的列i，其中包含观测类型“Type”。
读取SP3文件主要过程：
（1） 打开文件，然后读取文件，利用ReadSP3.m以矩阵的形式返回卫星的轨道信息、卫星数目、及文件头中的各种参数值。
（2） 读取精密轨道信息。利用for循环，得到卫星的编号矩阵以及其对应的轨道的位置及种差数据。
（3） 利用FormSP3.m根据指定的卫星编号以及起止时间读取SP3文件，并返回相关信息,以矩阵的形式输出。
涉及的主要函数及功能：
ReadSP3.m：读取SP3精密卫星轨道数据文件
调用方式：[SP3data,numsat,header]= ReadSP3(filename)
SP3data:轨道信息
Numsat：卫星数目
Header：文件头
FormSP3.m：读取SP3精密卫星轨道数据文件，并按卫星编号存为矩阵形式。
调用方式：[Time,SP3X,SP3Y,SP3Z,Clk,remark,new_prn]= FormSP3(filename,PRN,sep,lep)根据指定的卫星编号以及起止时间读取SP3文件filename，并返回相关信息。
四、实验体会
本实验主要是学习对N文件、SP3文件、O文件的读取，在实验过程中我学习到了很多数据的读取方法，这个实验之前我也遇到过好多数据读取方面的问题，学习了这个程序代码后，发现它读取数据的很多思路非常的巧妙。通过本次实验，各种文件格式数据的结构有
了深入的了解。
实验三卫星轨道计算
主要过程：
整个程序的实现分三个模块实现：时间转换模块；读广播星历文件(N文件)模块；计算卫星位置模块。
1、时间转换模块：
将O文件的某一历元对应的格里高利历时(实验中的时间：2001,9,4,9+40/60)转换为儒略时，然后将儒略时转换为GPS时得到对应的周数(WN)和周内秒数(TOW)。
主要涉及的函数及功能：
1）Julday.m：将格里高利历时转换为儒略时。
2）gps_time.m:将儒略时转换为GPS时的到对应的周数和周内秒数。
2、读星历文件模块：
在求解卫星位置时，第一需要利用O文件中每个历元的历元时刻t。在计算某时刻卫星的位置时，这里的某时刻便是O文件历元时刻t。第二需要利用读取的每个历元不同的卫星PRN号。根据PRN号和历元时刻t在广播星历N文件中搜索相同的卫星PRN号、合适的历元时刻，利用其对应的数据，计算卫星的位置。
先利用rinexe.m函数读取N文件将导航信息中的每颗卫星的21个参数以矩阵的格式存入文件eph.dat中，然后利用get_eph.m函数将文件eph.dat中的星历信息存入矩阵Eph中。
主要涉及的函数及功能：
1）rinexe.m：读取rinexe导航信息文件，将信息以矩阵的格式写入矩阵。矩阵中有21行，每一列存储一颗卫星的信息。
调用方式：rinexe(ephemerisfile, outputfile)，
ephemerisfile为文件名(N文件要放在rinexe函数所在的目录下)
outputfile为文本文件名，将读到N文件的数据存放在此文本文件中。
2）get_eph.m:星历文件中包含的星历信息被存入矩阵eph中。矩阵行数为21行，列数为有星历信息的卫星个数。
调用方式: eph = get_eph(ephemeridesfile)
Ephemeridesfile——文件名，将存放在此文本文件中的导航电文数据提取出来，为计算卫星的位置做准备。该函数返回值为eph星历矩阵(为21*7的矩阵)，它为Matlab内部格式。
3、计算卫星位置模块：
主要用到四个函数，分别为读O文件函数、eph星历矩阵的选择函数、修复函数(修复GPS时间超限或者溢出)、计算卫星的位置函数。
根据要寻找的历元所包含的卫星PRN号，在N文件中搜索对应的卫星PRN号、合适的历元时刻，利用其对应的数据，计算出卫星的位置。
利用find_eph.m函数从eph矩阵中找出相应卫星编号的卫星的参数所在的列数，接着利用satpos.m函数来该颗卫星对应的X、Y、Z坐标值，即在地心地固坐标系中的坐标。
主要涉及的函数及功能：
1）find_eph.m：从Eph矩阵中选中用于计算的列数。
调用方式：col_Eph(t) = find_eph(Eph,svs(t),sow)
Eph——星历矩阵，存放计算卫星位置所需要的N文件的数据。
svs(t)——svs中存放O文件中某个历元观测到的所有卫星PRN号。t为循环控制，svs(t)为svs中的某刻卫星的编号。
Sow——选定的历元的周内秒数值。
返回值——col_Eph矩阵中存放着从Eph矩阵中选中用于计算的列数。
2）X=satpos(tx_GPS,Eph(:,k))
tx_GPS——历元的格里高利历时转化为GPS时后的秒数sow（即所谓的归化后的时间。
Eph(:,k)——Eph星历矩阵中的对应的列的数据。
返回值——卫星在地心地固坐标系中的坐标
在本实验中的具体调用方式：
sat(1:3,t) = satpos(sow,Eph(:,col_Eph(t)))
其中Eph(:,col_Eph(t))为上一步中用col_Eph.m函数并找到的相应卫星在Eph矩阵中对应列的数据。
3）tk = check_t(t-toe)
t-toe——为儒略日
返回值为修复后的儒略日。
4）计算卫星在地固坐标系中的三维坐标的在程序中的具体实现如下： satp(1,1) = x1*cos(Omega)-y1*cos(i)*sin(Omega);
satp(2,1) = x1*sin(Omega)+y1*cos(i)*cos(Omega);
satp(3,1) = y1*sin(i);%计算卫星地心固定坐标系中的直角坐标 附：
1、在编计算每颗卫星位置的程序时，每颗卫星都需要存储导航电文中的21个参数，分别为：
svprn,af2,M0,roota,deltan,ecc,omega,cuc,cus,crc,i0,idot,cic,cis,Omega0,Omegadot,toe,af0,af1,toc。详细情况可以参照代码进行调试理解。
2、实验内容中要求的第3条“根据卫星位置计算公式编写主函数，同时调用时间转换模块、星历读取等的子函数来共同完成卫星位置的计算，最后输出结果”，具体的实现情况可以参照已经给出的代码。
实验四 伪距定位及DOP值
一、实验目的
1、掌握伪距离定位的基本原理并了解其公式推导基本流程
2、熟练掌握线性化的基本公式及内涵
3、掌握DOP计算的基本步骤并绘制时段内星座DOP的分布
二、实习内容
1、独立调试一段线性化程序，程序的输出为A矩阵，用于最小二乘求解与DOP值的计算；
2、构建A矩阵并提取对应DOP计算所需元素
3、空间直角坐标系转换至大地直角坐标系的R矩阵的编写
4、绘制指定接收机位置及指定时段内的星座DOP分布图
三、实验过程
1、主程序功能模块
（1） 调用函数mpgen计算多路径效应误差；
（2） 调用函数llh2xyz将用户坐标转换为WGS4坐标系统下的坐标；
（3） 调用函数loadgps加载gps数据；
（4） 调用函数olspos计算方向余弦矩阵H；
（5） 调用函数plot绘制DOP值图像
2、olspos函数功能模块
根据伪距观测值和卫星坐标计算用户坐标和方向余弦矩阵H；
该函数还利用循环进行迭代平差；
while ( (itertol) )
for N = 1:numvis
pr0 = norm(svxyzmat(N,:)-estusr(1:3));
y(N,1) = prvec(N)estusr(4);
end
H = hmat(svxyzmat,estusr(1:3));
beta = H\\y;
estusr=estusr+beta\';
iter=iter+1;
end
3、hmat函数模块
根据伪距方程线性化后的函数，求出三个该正数的系数，然后组建h矩阵;
(Xj-X)2(Yj-Y)2(Zj-Z)2
(Xj,Yj,Zj)为第j颗卫星在空中的位置；
(X,Y,Z)为接收机空间的位置。
线性化后得到：
Xj-X0Yj-Y0Zj-Z0 00dx0dy0dz
四、实验体会
通过本次实验的学习，我对伪距离定位的基本原理并及其公式推导基本流程有了更清楚的了解。同时也掌握了线性化的基本公式及内涵、精度因子计算的基本步骤，能够绘制了时段内星座DOP的分布图。
