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

卫星导航定位算法与程序设计实验报告.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文件读取，然后后面的计算代码较为复杂，需要对卫星轨道计算的公式了解全面，因此对这部分实验，是在老师和同学的帮助下完成的，直接在轨道计算的函数下加入了输出函数，将计算的轨道数据直接输出出来。

通过这部分实验是我加深了对卫星轨道计算公式的了解，以及以后面对复杂的公式应该如何应对。

。