卫星导航定位算法与程序设计实验报告.docx
《卫星导航定位算法与程序设计实验报告.docx》由会员分享，可在线阅读，更多相关《卫星导航定位算法与程序设计实验报告.docx（15页珍藏版）》请在冰豆网上搜索。
卫星导航定位算法与程序设计实验报告
2013级测绘工程专业
卫星导航定位算法与程序设计
实
验
报
告
实验名称：
卫星导航基本程序设计
班级：
学号：
姓名：
实验时间：
2016年6月28日~2016年6月30
中国矿业大学
实验一时空基准转换
一、实验目的
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”的头文件，建立类“CTimeTr”，创建变量“GPSTime”、“Time”、”JulDay”，并且申明函数“TIME2JUL”、“TIME2GTIME”等，用这些函数分别实现所需要的转换。
1.2具体的实现函数：
“TIME2JUL”函数：
doubleCTimeTr:
:
TIME2JUL（）//TIMETime,JULIANDAY&JulDay
{
doublem,y;
doubleD;
//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）;
return0;
}
“TIME2GTIME”：
voidCTimeTr:
:
TIME2GTIME（）
{
doubleJD;
longm,y;
intWN;
doubleWsecend;
//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;
whileabs（B1-B0）>=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*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;
end
L=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”函数：
functioneph=get_eph（ephemeridesfile）
fide=fopen（ephemeridesfile）;
[eph,count]=fread（fide,Inf,'double'）;
noeph=count/22;
eph=reshape（eph,22,noeph）;
“rinexe”函数（部分）：
functionrinexe（ephemerisfile,outputfile）
fide=fopen（ephemerisfile）;
head_lines=0;
while1
head_lines=head_lines+1;
line=fgetl（fide）;
answer=findstr（line,'ENDOFHEADER'）;
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”函数，将所需要的数据依次打印出来。
2.2、主要函数：
“anheader”函数：
function[Obs_types,ant_delta,ifound_types,approx_XYZ]=anheader（file）
fid=fopen（file,'rt'）;
eof=0;
ifound_types=0;
Obs_types=[];
ant_delta=[];
approx_XYZ=[];
while1%Gobblingtheheader
line=fgetl（fid）;
answer=findstr（line,'ENDOFHEADER'）;
if~isempty（answer）,break;end;
if（line==-1）,eof=1;break;end;
answer=findstr（line,'ANTENNA:
DELTAH/E/N'）;
if~isempty（answer）
fork=1:
3
[delta,line]=strtok（line）;
del=str2num（delta）;
ant_delta=[ant_deltadel];
end;
end
answer=findstr（line,'APPROXPOSITIONXYZ'）;
if~isempty（answer）
fork=1:
3
[app_XYZ,line]=strtok（line）;
del=str2num（app_XYZ）;
approx_XYZ=[approx_XYZdel];
end;
end
answer=findstr（line,'#/TYPESOFOBSERV'）;
if~isempty（answer）
[NObs,line]=strtok（line）;
NoObs=str2num（NObs）;
fork=1:
NoObs
[ot,line]=strtok（line）;
Obs_types=[Obs_typesot];
end;
ifound_types=1;
end;
end;
%fclose（fid）;
“grabdata”函数：
functionObs=grabdata（fid,NoSv,NoObs）
%GRABDATAPositionedinaRINEXfileataselectedepoch
%readsobservationsofNoSvsatellites
globallin
Obs=zeros（NoSv,NoObs）;
ifNoObs<=5%ThiswilltypicalbeTurboSIIdata
foru=1:
NoSv
lin=fgetl（fid）;
fork=1:
NoObs
Obs（u,k）=str2num（lin（2+16*（k-1）:
16*k-2））;
end
end
else%ThiswilltypicalbeZ12data
Obs=Obs（:
[12345]）;%Wecancelthelasttwocolumns6and7
NoObs=5;
foru=1:
NoSv
lin=fgetl（fid）;
lin_doppler=fgetl（fid）;
fork=1:
NoObs%%-1
ifisempty（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））;
end
end
end
2.3实验结果
四、实验感想
这部分实验是我在之前做的，之前自己有看过“gps_easy”有关的代码，看过相关的“N文件”“O文件”读写函数，并且学会了如何调用这些函数，对里面的输出量有了一点的了解，所以我自己的主要工作就是运用了“fprintf”函数，将读取到matlab中的矩阵写入TXT文档中，这部分工作量不是很大，但较有意义。
实验三卫星轨道计算
一、实验目的
1、进一步熟悉N文件的读入
2、掌握开普勒参数计算卫星轨道的过程
3、编程实现采用广播星历计算卫星轨道
4、掌握MATLAB函数调用步骤
二、实验内容
本实验内容包括：
1、调试时间转换函数，熟悉内容，备主函数调用
2、调试广播星历导航文件的读入程序，备主函数调用
3、根据卫星位置计算公式编写主函数，同时调用时间转换、星历读取等的子函数来共同完成卫星位置的计算，最后输出计算结果
4、理清程序各模块的功能结构
三、实验过程
1、实验思路：
在老师提供的“SPP”文件中，直接利用卫星位置计算函数，进行卫星位置的计算，将利用“Gps.cpp”文件中的”GetGpsPosition”函数，利用其中的迭代求解卫星位置部分，用“cout”直接将卫星迭代后的位置直接输出，因为星历文件中有较多的星历，所以利用循环语句，将求解出来的卫星位置依次输出出来。
2、主要函数
boolCGps:
:
GetGpsPosition（）
{
GPSTIMEts;
GPSTIMEtr;
GPSTIMEts0;
GPSTIMEoGTime;
GPSTIMEnGTOC;
vectorSendSignPosition;
GpsPosGpsPTemp;
GpsSendPositionSdSignPoTemp;
intnTheFitPoint=0;
PositionpTemp;
GpsPosition.clear（）;
if（oData.size（）==0）
{
strErr=（"PRN=%d没有对应星历"）;
returnfalse;
}
if（nData.size（）==0）
{
returnfalse;
}
cout<for（inti=0;i{
double*VX=newdouble[4];//
//
for（intj=0;j{
//
TimeToGpsTime（oData[i].ObserveInfo.ObserveTime,oGTime）;
tr=oGTime;//
//
FindTheBestFitTime（oData[i].ObserveInfo.ObserveTime,
oData[i].oObserveData[j].PRN,nTheFitPoint）;
if（nTheFitPoint==-1）
{
//stringjuge;
//juge.（"PRN=%d没?
有®D对?
应®|星?
历¤¨²",oData[i].oObserveData[j].PRN）;
//AfxMessageBox（juge）;
continue;
}
ts.lWeek=tr.lWeek;
ts0.lWeek=tr.lWeek;
//
ts.lSecond=tr.lSecond-0.075;
do
{
ts0.lSecond=ts.lSecond;
SatellitePosition（nData[nTheFitPoint],ts0,pTemp）;
ts.lSecond=tr.lSecond-sqrt（pow（（pTemp.XX-x0）,2）+pow（（pTemp.YY-y0）,2）
+pow（（pTemp.ZZ-z0）,2））/c;
}while（fabs（ts0.lSecond-ts.lSecond）>1e-007）;//
//
earthrot（tr.lSecond-ts.lSecond,pTemp.XX,pTemp.YY,pTemp.ZZ）;
TimeToGpsTime（nData[nTheFitPoint].TOC,nGTOC）;
//
SdSignPoTemp.dt=nData[nTheFitPoint].dClkBias+
nData[nTheFitPoint].dClkDrift*（ts.lSecond-nGTOC.lSecond）+
nData[nTheFitPoint].dClkDriftRate*
pow（（ts.lSecond-nGTOC.lSecond）,2）;
//
SdSignPoTemp.PRN=oData[i].oObserveData[j].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=oData[i].oObserveData[j].dC1;
if（i<10）
cout<<<""<//
SendSignPosition.push_back（SdSignPoTemp）;
}
}
returntrue;
}
3、实验结果
四、实验感想
计算卫星位置的这部分“C++”代码比较综合，读代码的难度较大，它综合了前面的“C++”O文件读取，N文件读取，然后后面的计算代码较为复杂，需要对卫星轨道计算的公式了解全面，因此对这部分实验，是在老师和同学的帮助下完成的，直接在轨道计算的函数下加入了输出函数，将计算的轨道数据直接输出出来。
通过这部分实验是我加深了对卫星轨道计算公式的了解，以及以后面对复杂的公式应该如何应对。
。
升级会员 