卫星导航定位算法与程序设计实验报告.docx
《卫星导航定位算法与程序设计实验报告.docx》由会员分享,可在线阅读,更多相关《卫星导航定位算法与程序设计实验报告.docx(16页珍藏版)》请在冰豆网上搜索。
![卫星导航定位算法与程序设计实验报告.docx](https://file1.bdocx.com/fileroot1/2023-1/23/6812219d-d56f-43a0-95dc-e5883d8463ef/6812219d-d56f-43a0-95dc-e5883d8463ef1.gif)
卫星导航定位算法与程序设计实验报告
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、每人独立完成实验规定的内容
二、实验内容
本实验内容包括:
内容一:
编程实现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++”代码比较综合,读代码的难度较大,它综合了前