卫星导航定位算法与程序设计实验报告.docx
- 文档编号:29408384
- 上传时间:2023-07-23
- 格式:DOCX
- 页数:19
- 大小:124.94KB
卫星导航定位算法与程序设计实验报告.docx
《卫星导航定位算法与程序设计实验报告.docx》由会员分享,可在线阅读,更多相关《卫星导航定位算法与程序设计实验报告.docx(19页珍藏版)》请在冰豆网上搜索。
卫星导航定位算法与程序设计实验报告
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”:
void CTimeTr:
:
TIME2GTIME()
{
doubleJD;ﻩ
ﻩ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;
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 %Gobbling theheader
line= fgetl(fid);
answer=findstr(line,’END OFHEADER');
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_types ot];
end;
ifound_types=1;
end;
end;
%fclose(fid);
“grabdata”函数:
functionObs=grabdata(fid, NoSv,NoObs)
%GRABDATAPositionedina RINEX a selectedepoch
% reads observations of NoSv satellites
globallin
Obs=zeros(NoSv,NoObs);
ifNoObs〈=5 %This willtypicalbe TurboSIIdata
foru= 1:
NoSv
lin=fgetl(fid);
fork =1:
NoObs
Obs(u,k) = str2num(lin(2+16*(k—1):
16*k-2));
end
end
else %ThiswilltypicalbeZ12 data
Obs=Obs(:
[1 2345]);%Wecancelthelasttwocolumns 6and 7
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;
ﻩGPSTIME oGTime;
ﻩGPSTIMEnGTOC;
ﻩvector ﻩGpsPosGpsPTemp; GpsSendPositionSdSignPoTemp; intnTheFitPoint=0; PositionpTemp;ﻩ ﻩGpsPosition.clear();ﻩ ﻩif (oData、size()==0) ﻩ{ strErr= (”PRN=%d没有对应星历”); ﻩreturnfalse; ﻩ} ﻩif (nData。 size()==0) { returnfalse; ﻩ} ﻩcout<〈fixed; ﻩfor(int i=0;i〈oData、size();i++) ﻩ{ ﻩdouble*VX=newdouble[4];// ﻩ// ﻩfor(intj=0;j〈oData[i]、ObserveInfo。 ObserveSatSum;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< PRN〈<"”< ZZ ﻩ<〈” "< SendTime。 lWeek〈〈””<
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 卫星 导航 定位 算法 程序设计 实验 报告