卫星导航定位实验报告

China University of Mining and Technology 《卫星导航定位算法与程序设计》

实验报告

学号： 07122825

姓名：王亚亚

班级：测绘12—1

指导老师：王潜心/张秋昭/刘志平

中国矿业大学环境与测绘学院

2015-07-01

实验一编程实现读取下载的星历

一、实验要求：

读取RINEX N 文件，将所有星历放到一个列表（数组）中。并输出和自己学号相关的卫星编号的星历文件信息。读取RINEX O文件，并输出指定时刻的观测信息。

二、实验步骤：

1、下载2014年的广播星历文件和观测值文件，下载地址如下：

ftp://https://www.360docs.net/doc/2c10762235.html,/gps/data/daily/2014/

2、要求每一位同学按照与自己学号后三位一致的年积日的数据文件和星历文件，站点的选择必须选择与姓氏首字母相同的站点的数据，以王小康同学为例，学号：07123077，需下载077那天的数据。有些同学的学号365<后三位

<730，则取学号后三位-365，以姜平同学为例：学号10124455，下载455-

365=90 天的数据，有些同学的学号730<后三位<=999，则取学号后三位-730，以万伟同学为例：学号：07122854，则下载854-730 = 124天的数据。可以选择wnhu0124.14n wnhu0124.14o 根据上述要求我下载了2014年第95天的数据，选择其中的wsrt0950.14n和wsrt0950.14o星历文件。指定时刻（学号后五位对应在年积日对应的秒最相近时刻）的观测值信息如张良09123881，后五位23881，取23881-3600*6= 2281秒，6点38分01秒，最近的历元应该是6点38分00秒的数据。根据计算与我最接近的观测时刻为2014年4月5日6点20分30.00秒。

3、编程思路：

利用rinex函数读取星历文件中第14颗卫星的星历数据并输出显示。对数据执行762次循环找到对应的2014年4月5日6点20分30.00秒，并输出观测值。

4、程序运行结果：

三、编程主要代码：

1、读写N 文件主函数代码： clc;%清除屏幕之前打过的代码

clear all ;%从内存中释放所有的内存变量和数组以及所有用户自定义的菜单栏, 菜单, 和窗口的定义。

close all ;%关闭当前及所有工作区中所有打开的数据库, 表, 和索引，并选择 1 号工作区 [filename,filepath]=uigetfile('*.14n','选择计算的输入文件');%打开窗口选取计算的输入文件 file=[filepath filename];%获得文件名 rinexe(file,'eph.dat');%函数rinexe 读取一个标准格式化导航信息文件并重新格式化数据为21行列数等于卫星数的矩阵。

Eph = get_eph('eph.dat');%将存放在‘eph.dat ’中的导航电文数据提取出来，为计算卫星位置做准备，返回值-eph 星历矩阵，将读到的数据存入内存（注意：eph 星历矩阵，每颗卫星存放21个参数） Eph

[filename,filepath]=uiputfile('王亚亚-作业1读写N 导航文件.txt','选择一个路径对结果文件保存');%选择路径保存输出的结果文件 file=[filepath filename];%获取文件名 fid = fopen(file,'wt+');%文件有两种格式：二进制文件(b)和文本文件(t)。在Windows 下，打开文件的默认是二进制格式，如果要以文本方式打开，则必须在打方式中加上字符 't'。 fprintf(fid,'%c',' 班级 姓名 学号 作业序号 主要内容 ');

fprintf(fid,'%c\n',' ');

fprintf(fid,'%c',' 测绘12-1 王亚亚

07122825 上机实验作业一 读写N 导航文件 '); fprintf(fid,'%c\n',' '); for k = 1: size(Eph,1)

p=num2str(Eph(k,44));%输出第14颗卫星的星历数据

fprintf(fid,'%c',p);%以一个字符输出 if k == 1;

fprintf(fid,'%c',' %GPS 卫星PRN 编号'); end if k == 2;

fprintf(fid,'%c',' %GPS 卫星钟的漂移速度'); end

if k == 3; fprintf(fid,'%c',' %参考时刻的升交点赤经M0'); end

if k == 4; fprintf(fid,'%c',' %roota sqrt(A)( m1 / 2(根号下) ) 轨道长半轴的平方根'); end

if k == 5; fprintf(fid,'%c',' %deltan Δn(rad /s) 平均运动修正量'); end

if k == 6; fprintf(fid,'%c',' %ecc e 轨道偏心率'); end if k == 7;

fprintf(fid,'%c',' %omega ω(rad) 近地点的角距'); end

if k == 8; fprintf(fid,'%c',' %cuc Cuc(rad) 纬度幅角的余弦调和项改正的振幅'); end

if k == 9; fprintf(fid,'%c',' %cus Cus(radians) 纬度幅角的正弦调和项改正的振幅'); end

if k == 10; fprintf(fid,'%c',' %crc Crc(m) 轨道半径的余弦调和项改正的振幅');

end

if k == 11; fprintf(fid,'%c',' %crs Crs（m）轨道半径的正弦调和项改正的振幅'); end

if k == 12; fprintf(fid,'%c',' %i0 i0(rad) 轨道倾角'); end

if k == 13; fprintf(fid,'%c',' %idot i(rad/s)(IDOT) 轨道倾角变化率'); end

if k == 14; fprintf(fid,'%c',' %cic Cic(rad) 轨道倾角的余弦调和项改正的振幅'); end

if k == 15; fprintf(fid,'%c',' %cis

Cis(rad) 轨道倾角的正弦调和项改正的振幅'); end

if k == 16;

fprintf(fid,'%c',' %Omega0

Ω(rad)(OMEGA) 参考时刻的升交点的赤经'); end

if k == 17;

fprintf(fid,'%c',' %Omegadot Ω

(rad/s)(OMEGA DOT) 升交点赤经的变化率'); end

if k == 18;

fprintf(fid,'%c',' %toe TOE星历的参考时刻(GPS周内的秒数) 星历表参考历元'); end

if k == 19; fprintf(fid,'%c',' %af0 卫星钟的偏差（s）卫星钟偏差'); end

if k == 20; fprintf(fid,'%c',' %af1 卫星钟的漂移（s/s）'); end

if k == 21;

fprintf(fid,'%c',' %toe 存储TOE星历的参考时刻(GPS周内的秒数) 电文发送时刻'); end

fprintf(fid,'%c\n',' ');

end

fclose(fid);

2、读写O文件主函数代码：

clc;%清除屏幕之前打过的代码

clear all;%从内存中释放所有的内存变量和数组以及所有用户自定义的菜单栏, 菜单, 和窗口的定义。

close all;%关闭当前及所有工作区中所有打开的数据库, 表, 和索引，并选择1 号工作区% We identify the observation file and open it

我们确定观测文件并打开它

[filename,filepath]=uigetfile('*.14O','选择计算的输入的O文件');

file=[filepath filename];%文件名

fid1 = fopen(file,'rt');%r：读出观测文件ofile1，t：以文本的格式（因为fopen函数默认打开方式是二进制）

[Obs_types1, ant_delta1, ifound_types1, eof11] = anheader(file);%打开观测文件ofile1，从头文件中截取观测类型Obs_types1、天线参数的信息ant_delta1，以及标明是否寻找到观测类型行ifound_types1、是否到文件尾的判定eof11 NoObs_types1 = size(Obs_types1,2)/2;%返回观测类型Obs_types1列数(1:行2：列)的

Obs_types1矩阵是1行14列：[C1P1P2L1L2D1D2] size（Obs_types1,2）值是14除以2得观测值类型个数是7种

Pos = [];

% There are 20 epochs of data in ofile1

SITE247J.01N文件中有20个历元的数据

for q = 1:761 %循环762次找到对应的2014年4月5日6点20分30.0000000秒

[time1, dt1, sats1,eof1,datee] =

fepoch_0(fid1);

% [time1, dt1, sats1, eof1] =

fepoch_0(fid1);%通过函数fepoch_0调用观测文件SITE247J.01O，返回四个参数：time1：周内秒、dt1：接收机时钟的偏差( 单位为s, 为可选项)、sats1：卫星的PRN号矩阵、

eof1：返回是否到文件尾的标示符（1：结束0：未结束）

NoSv1 = size(sats1,1);%通过查看行向量sats1的长度，获取卫星个数

% We pick the observed P2 pseudoranges 我们获取观测的P2码伪距其中伪随机噪声码（PN）即测距码主要有精测距码（P码）和粗测距码（C/A码）两种。其中P码的码率

为10.23MHz、C/A码的码率为1.023MHz，长度为1023

obs1 = grabdata(fid1, NoSv1,

NoObs_types1); %输入观测文件名、观测卫星个数、观测卫星类型数。输出Obs矩阵，行：卫星个数、列：观测类型：C1 P1 P2 L1 L2 D1 D2

i = fobs_typ(Obs_types1,'P2');%在观测类型：C1 P1 P2 L1 L2 D1 D2查找出‘P2’所在位置，并赋值给i

2.235646758400000e+07

end

[filename,filepath]=uiputfile('王亚亚-作业1读写O观测文件.txt','选择一个路径对结果文件保存');%选择路径保存输出的结果文件

file=[filepath filename];

fid = fopen(file,'wt+');%文件有两种格式：二进制文件(b)和文本文件(t)。在Windows下，打开文件的默认是二进制格式，如果要以文本方式打开，则必须在打方式中加上字符't'。

%fid = fopen('作业1.txt','wt+');%文件有两种格式：二进制文件(b)和文本文件(t)。在Windows 下，打开文件的默认是二进制格式，如果要以文本方式打开，则必须在打方式中加上字符't'。

fprintf(fid,'%c',' 班级姓名

学号作业序号主要内容');

fprintf(fid,'%c\n',' ');

fprintf(fid,'%c',' 测绘12-1 王亚亚07122825 上机实验作业一读写观测O文件');fprintf(fid,'%c\n',' ');

%fprintf(fid,'%8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n',datee');

for i=1:1%输出datee矩阵即时间

for j=1:6

fprintf(fid,'%d\t',datee(i,j));

if rem(i,1)==0 && rem(j,6)==0

fprintf(fid,'%c','%时间');

%fprintf(fid,'\r\n');

end

end

end

fprintf(fid,'%c\n',' ');

fprintf(fid,'%c',' L1 L2

C1 P1 P2 ');

fprintf(fid,'%c\n',' ');

for k=1:size(sats1,1);%卫星个数

f=sats1(k,1);

for i=1:1%输出obs1矩阵

for j=1:5

fprintf(fid,'%d\t',obs1(i,j));

if rem(i,1)==0 && rem(j,5)==0

fprintf(fid,'%s%d%s ','%读取的是',f','号卫星');

fprintf(fid,'\r\n');

end

end

end

end

fprintf(fid,'%c\n',' ');

fclose(fid);

3、主要函数代码及其注释：

rinex函数代码：

function rinexe(ephemerisfile, outputfile)

fide = fopen(ephemerisfile); %fopen（）函数功能是打开一个文件，在这里就是打开星历表文件

head_lines = 0; %头行赋值为0

while 1 % We skip header 我们跳过头行

那么while 1 其中1代表一个常量表达式，他永远不会等于0。所以，循环会一直执行下去。除非你设置break等类似的跳出循环语句循环才会中止。该循环用于跳过整个头文件head_lines = head_lines+1;%当跳出循环体的时候，head_lines存储了头文件的行数为多少

line = fgetl(fide); %fgetl（）函数功能：在MATLAB中，该函数用于从文件中读取一行数据，并丢弃行末的换行符；读取一行后停

留在该行末端，line保存的是字符串

answer = findstr(line,'END OF

HEADER'); %findstr（）在文件中寻找字符串END OF HEADER 注释RINEX 格式文件的文件头用于存放与整个文件有关的全局性信息, 位于每个文件的最前部, 其最后一个记录为“END OF HEADER”。

if ~isempty(answer), break; end; %如果answer是非空的，则跳出当前的循环体end;

head_lines

noeph = -1;

while 1

noeph = noeph+1;%循环结束时存储了星历表的行数

line = fgetl(fide);

if line == -1, break; end%line是一个文本字符串，除非该行只包含结束的文件标记。在这种情况下，tline是数字值-1。返回由文件标识符fid指示的文件的下一行。如果fgetl遇到文件结束指示符，则返回-1

end;

noeph = noeph/8 %对观测数据进行分组，8行数据为一组每颗卫星一个参考时刻的数据占8 行, 第1 行为卫星的PRN 号和该卫星时钟的参考时刻及其改正模型参数, 第2～8 行为该卫星的广播轨道数据。

frewind(fide); %位置指针移至文件fide首部for i = 1:head_lines, line = fgetl(fide); end; %提取头文件的那几行

% Set aside memory for the input %设置另外一个存储体用于输入

svprn = zeros(1,noeph); %创建一行N （卫星颗数）列的向量空间存储卫星导航数据

weekno = zeros(1,noeph);

t0c = zeros(1,noeph);

tgd = zeros(1,noeph);

aodc = zeros(1,noeph);

toe = zeros(1,noeph);

af2 = zeros(1,noeph);

af1 = zeros(1,noeph);

af0 = zeros(1,noeph);

aode = zeros(1,noeph);deltan = zeros(1,noeph);

M0 = zeros(1,noeph);

ecc = zeros(1,noeph);

roota = zeros(1,noeph);

toe = zeros(1,noeph);

cic = zeros(1,noeph);

crc = zeros(1,noeph);

cis = zeros(1,noeph);

crs = zeros(1,noeph);

cuc = zeros(1,noeph);

cus = zeros(1,noeph);

Omega0 = zeros(1,noeph);

omega = zeros(1,noeph);

i0 = zeros(1,noeph);

Omegadot = zeros(1,noeph);

idot = zeros(1,noeph);

accuracy = zeros(1,noeph);

health = zeros(1,noeph);

fit = zeros(1,noeph);

for i = 1:noeph

line = fgetl(fide); %接着上面的提取卫星星历表第一行数据记录，前面几行是头文件。第一行信息：《《《《PRN 号/ 历元/ 卫星钟》》》》

svprn(i) = str2num(line(1:2));%把字符串前2个字符转换数值（PRN号），可以使用

str2double代替str2num。在matlab中，无论是内建函数还是工具箱函数，2很常见，这可能是因为2英文two和to发音相同。而2写起来也比较简单。所以很多转换类函数都用2来命名而非to。第1 行为卫星的PRN 号和该卫星时钟的参考时刻及其改正模型参数

year = line(3:6);%把字符第3至6的字符转化为年（数值型）

month = line(7:9);%把字符第7至9的字符转化为月（数值型）

day = line(10:12);%把字符第10至12的字符转化为日（数值型）

hour = line(13:15);%把字符第13至15的字符转化为时（数值型）

minute = line(16:18);%把字符第16至18的字符转化为分（数值型）

second = line(19:22);%把字符第19至22的字符转化为秒（数值型）

af0(i) = str2num(line(23:41));%卫星钟的偏差（s）就是第一行的第二段

af1(i) = str2num(line(42:60));%卫星钟的漂移（s/s）就是第一行的第三段

af2(i) = str2num(line(61:79));%卫星钟的漂移速度(s/s2平方)就是第一行的第四段

line = fgetl(fide); %读取的是卫星星历表第二行数据记录。第二行信息：《《《《广播轨道—1》》》》

IODE = line(4:22);%IODE( Issue of Data, Ephemeris /数据、星历发布时间)就是第二行的第一块

crs(i) = str2num(line(23:41));%Crs（m）第二行的第二段

deltan(i) = str2num(line(42:60));%Δn(rad /s)第二行的第三段

M0(i) = str2num(line(61:79));%M0下标( rad)第二行的第四段

line = fgetl(fide); %读取的是卫星星历表第三行数据记录。第三行信息：《《《《广播轨道—2》》》》

cuc(i) = str2num(line(4:22)); %Cuc(rad)

ecc(i) = str2num(line(23:41));%e轨道偏心率

cus(i) = str2num(line(42:60));%Cus(radians)

roota(i) = str2num(line(61:79));%sqrt(A)( m1 / 2(根号下) )

line=fgetl(fide); %读取的是卫星星历表第四行数据记录。第四行信息：《《《《广播轨道—3》》》》

toe(i) = str2num(line(4:22));%TOE星历的参考时刻(GPS周内的秒数)

cic(i) = str2num(line(23:41));%Cic(rad)

Omega0(i) =

str2num(line(42:60));%Ω(rad)(OMEGA)

cis(i) = str2num(line(61:79));%Cis(rad)

line = fgetl(fide); %读取的是卫星星历表第五行数据记录。第五行信息：《《《《广播轨道—4》》》》

i0(i) = str2num(line(4:22));%i0(rad)

crc(i) = str2num(line(23:41));%Crc(m)

omega(i) = str2num(line(42:60));%ω(rad)

Omegadot(i) =

str2num(line(61:79));%Ω(rad/s)(OMEGA DOT) line = fgetl(fide); %读取的是卫星星历表第六行数据记录。第六行信息：《《《《广播轨道—5》》》》

idot(i) = str2num(line(4:22));%i(rad/s)(IDOT) codes = str2num(line(23:41));%L2上的码

weekno = str2num(line(42:60));%PS 周数( 与TOE 一同表示时间)。为连续计数,不是1024的余数

L2flag = str2num(line(61:79));%L2 P码数据标记

line = fgetl(fide); %读取的是卫星星历表第七行数据记录。第七行信息：《《《《广播轨道—6》》》》

svaccur = str2num(line(4:22));%卫星精度( m) svhealth = str2num(line(23:41));%卫星健康状态(第1子帧第3字第17～22 位)

tgd(i) = str2num(line(42:60));%TGD(sec)

iodc = line(61:79);%IODC钟的数据龄期

line = fgetl(fide); %读取的是卫星星历表第八行数据记录。第八行信息：《《《《广播轨道—7》》》》

tom(i) = str2num(line(4:22));%电文发送时刻(单位为GPS周的秒, 通过交接字(HOW)中的Z 计数得出)

spare = line(23:41);%拟合区间(h),如未知则为零

% spare = line(42:60);%备用

% spare = line(61:79);%备用

end

status = fclose(fide) %关闭文件fide

% Description of variable eph. 星历矩阵变量的描述（每颗卫星有21个变量，列为卫星颗数？21颗0 0？）

eph(1,:) = svprn; %星历表矩阵第一行存储每颗卫星的PRN号（SVs号段）

eph(2,:) = af2;%第二行存储卫星钟的漂移速度

eph(3,:) = M0;

eph(4,:) = roota;

eph(5,:) = deltan;

eph(6,:) = ecc;

eph(7,:) = omega;

eph(8,:) = cuc;

eph(9,:) = cus;

eph(10,:) = crc;

eph(11,:) = crs;

eph(12,:) = i0;

eph(13,:) = idot;

eph(14,:) = cic;

eph(15,:) = cis;

eph(16,:) = Omega0; eph(17,:) = Omegadot; eph(18,:) = toe;

eph(19,:) = af0;

eph(20,:) = af1;eph(21,:) = toe;%第21行存储TOE星历的参考时刻(GPS周内的秒数)

fidu = fopen(outputfile,'w');% w 打开只写文件，若文件存在则文件长度清为0，即该文件内容会消失。若文件不存在则立刻创建文件count = fwrite(fidu,[eph],'double');%将星历矩阵中的元素写到指定文件fidu中去，精度类型为双精度double

fclose all%关闭所有文件

四、实验总结：

通过本次实验我们更加清晰的了解到了卫星星历文件的数据信息及其含义。对卫星的导航文件的格式信息也有了深刻的了解。在编程中遇到的最大问题是，习惯了C++编程的方法，此次用Matlab进行编程经常遇到各种错误信息，尤其是格式转换问题，经常要用到‘movefile’函数对数据进行格式转换。经过此次实验深刻体会到批处理和利用程序代替人工执行工作的效率。为以后从事和GPS 有关的工作奠定了一定基础。

实验二输出GPS卫星坐标值

一、实验要求：

利用下载的对应星历文件计算指定时刻（学号后五位对应在年积日对应的秒）的所有卫星的位置，如张良09123881，后五位23881，取23881-3600*6= 2281秒，6点38分01秒。根据我的学号计算后为读取2014年4月5日6点20分30.00秒的星历文件。并将此时段的卫星坐标信息输出。

二、实验步骤：

1、数据准备，采用实验一下载2014年的广播星历文件和观测值文件。

2、找到GPS卫星号并把卫星坐标输出。

3、编程思路：

通过卫星观测文件找到卫星编号，并将卫星的坐标值匹配输出。

4、程序运行结果：

三、编程主要代码：

1、主函数代码：

clc;%清除屏幕之前打过的代码

clear all ;%从内存中释放所有的内存变量和数组以及所有用户自定义的菜单栏, 菜单, 和窗口的定义。

jd = julday(2014,4,5,6+40/60);

[week,sow] = gps_time(jd);

[filename,filepath]=uigetfile('*.14n','选择计算的输入文件');%打开窗口选取计算的输入文件file=[filepath filename];%获得文件名

rinexe(file,'eph.dat');%函数rinexe读取一个标注格式化导航信息文件并重新格式化数据为21行列数等于卫星数的矩阵。

Eph = get_eph('eph.dat');

svs = [1 4 10 11 12 13 17 20 23 25 31 32];%观测文件卫星

for t = 1:length(svs)

col_Eph(t) = find_eph(Eph,svs(t),sow);

sat(1:3,t) = satpos(sow,Eph(:,col_Eph(t))); end

sat % position of svs in ECI system %%%%%%%%%%%%%%%%%%%%% end easy2.m %%%%%%%%%%%%%%% [filename,filepath]=uiputfile('王亚亚-作业2输出GPS卫星XYZ位置.txt','选择一个路径对结果文件保存');%选择路径保存输出的结果文件

file=[filepath filename];

fid = fopen(file,'wt+');%文件有两种格式：二进制文件(b)和文本文件(t)。在Windows下，打开文件的默认是二进制格式，如果要以文本方式打开，则必须在打方式中加上字符't'。fprintf(fid,'%c',' 班级姓名

学号作业序号主要内容');

fprintf(fid,'%c\n',' ');

fprintf(fid,'%c',' 测绘12-1 王亚亚07122825 上机实验作业二输出GPS 卫星XYZ位置');

fprintf(fid,'%c\n',' ');

for k = 1: size(sat,2)

p=num2str(sat(1,k));

fprintf(fid,'%c',p);

f=svs(1,k);

fprintf(fid,' %s%d%s ','%GPS',f','号卫星X坐标');

fprintf(fid,'%c\n',' ');

p=num2str(sat(2,k));

fprintf(fid,'%c',p);

fprintf(fid,' %s%d%s ','%GPS',f','号卫星Y坐标');

fprintf(fid,'%c\n',' ');

p=num2str(sat(3,k));

fprintf(fid,'%c',p);

fprintf(fid,' %s%d%s ','%GPS',f','号卫星Z坐标');

fprintf(fid,'%c\n',' ');

end

fclose(fid);

四、实验总结：

通过本次实验我们再一次认识到GPS卫星星历文件的读取，同时加深了我们对MATLAB语法的认识和简单的编程，也认识到大部分语言的逻辑思路上具有共通性。很庆幸的是MATLAB文件读取语法保留了C语言的文件读取格式，这样让编程的难度又降低了。同时也使我们清晰意识到以后编程的重要性。

9

实验三精密单点定位

一、实验要求：

编程实现GPS单点定位，要求加入电离层、对流层改正。

二、实验步骤：

1、数据准备，采用实验一下载2014年的广播星历文件和观测值文件。

2、读取2014年4月5日6点20分30.00秒到2014年4月5日6点30分30.00秒20次历元的星历文件数据，并计算wsrt0950观测站位置。

3、编程思路：

利用双频改正模型进行电离层延迟改正，利用Modified hopfield折射改正模式对对流层进行改正。

4、程序运行结果：

三、编程主要代码：

1、主函数代码：

clc ;%清除屏幕之前打过的代码

clear all;%从内存中释放所有的内存变量和数组以及所有用户自定义的菜单栏, 菜单, 和窗口的定义。

10

[filename,filepath]=uigetfile('*.14n','选择计算的N 文件');%打开窗口选取计算的输入文件

file=[filepath filename];%获得文件名

rinexe(file,'eph.dat');%函数rinexe读取一个标注格式化导航信息文件并重新格式化数据为21行列数等于卫星数的矩阵。

Eph = get_eph('eph.dat');%获得卫星星历表（Matlab的矩阵格式）

% We identify the observation file and open it

我们确定观测文件并打开它

[filename,filepath]=uigetfile('*.14o','选择计算的O 文件');%打开窗口选取计算的输入文件

file1=[filepath filename];%获得文件名

fid1 = fopen(file1,'rt');%r：读出观测文件

ofile1，t：以文本的格式（因为fopen函数默认打开方式是二进制）

[Obs_types1, ant_delta1, ifound_types1, eof11] = anheader(file1);%打开观测文件ofile1，从头文件中截取观测类型Obs_types1、天线参数的信息ant_delta1，以及标明是否寻找到观测类型行ifound_types1、是否到文件尾的判定eof11 NoObs_types1 = size(Obs_types1,2)/2;%返回观测类型Obs_types1列数(1:行2：列)的

Obs_types1矩阵是1行14列：

[C1P1P2L1L2D1D2] size（Obs_types1,2）值是14除以2得观测值类型个数是7种

Pos = [];%

for q = 1:762 %循环762次找到对应的2014年4月5日6点20分30.0000000秒

[time1, dt1, sats1,eof1,datee1] =

fepoch_0(fid1);

end

for q = 762:781%循环762次找到对应的2014年4月5日6点20分30.0000000秒

[time1, dt1, sats1, eof1,datee] =

fepoch_0(fid1);%通过函数fepoch_0调用观测文件SITE247J.01O，返回四个参数：time1：周内秒、dt1：接收机时钟的偏差( 单位为s, 为可选项)、sats1：卫星的PRN号矩阵、

eof1：返回是否到文件尾的标示符（1：结束0：未结束）

NoSv1 = size(sats1,1);%通过查看行向量sats1的长度，获取卫星个数

% We pick the observed P2 pseudoranges

我们获取观测的P2码伪距其中伪随机噪声码（PN）即测距码主要有精测距码（P码）和粗测距码（C/A码）两种。其中P码的码率为10.23MHz、C/A码的码率为1.023MHz，长度为1023

obs1 = grabdata(fid1, NoSv1,

NoObs_types1); %输入观测文件名、观测卫星个数、观测卫星类型数。输出Obs矩阵，行：卫星个数、列：观测类型：C1 P1 P2 L1 L2 D1 D2

i = fobs_typ(Obs_types1,'C1');%在观测类型：C1 P1 P2 L1 L2 D1

D2查找出‘P2’所在位置，并赋值给i

pos = recpo_ls(obs1(:,i),sats1,time1,Eph);%输入P2观测值（第i列）、卫星数、观测时的周内秒、星历矩阵

Pos = [Pos pos];%上面函数每计算一个pos （4行1列的pos存入Pos中，共循环了20次

q=1：20）一个q（历元数）列4行的矩阵，行：pos 为1、2、3：接收机位置4：接收机钟差El为5：高度角GDOP为6：GDOP basic_obs为7、8、9：改正后的坐标X 10：改正后的站星距离

end

me = mean(Pos,2); %对于矩阵mean(A,2)就是返回一个包含（每一行的平均值）的列向量fprintf('\nMean Position as Computed From 20 Epochs:') %由20个历元计算的平均位置fprintf('\n\nX: %12.3f Y: %12.3f Z: %12.3f T: %12.3f',me(1,1), me(2,1), me(3,1),me(4,1)) %输出X、Y、Z的平均值%12.3f这个是格式要求

figure(2);

plot((Pos(1:4,:)-Pos(1:4,1)*ones(1,q-

761))','linewidth',2) %ones（1，q）是生成一个1行q=20列的每个元素都为1的矩阵，则将矩阵所有坐标都减去了第一历元的卫星坐标，画图进行比较（循环完q=20）位置随时间变化图画出来

title('Positions Over Time','fontsize',16) %位置随时间的变化图% 画出其它历元与第一个历元xyz方向的差值

legend('X','Y','Z') %注释字符：'X','Y','Z'到图中

11

去

xlabel('Epochs [1 s interval]','fontsize',16) %x轴表示历元：一秒一间隔

ylabel('Variation in Coordinates, Relative to the First Epoch [m]','fontsize',16) %y坐标系相对于第一历元的变化

set(gca,'fontsize',16) %设置网格的显示格式legend %%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%

[filename,filepath]=uiputfile('王亚亚-作业3精密单点定位.txt','选择一个路径对结果文件保存');%选择路径保存输出的结果文件

file=[filepath filename];

fid = fopen(file,'wt+');%文件有两种格式：二进制文件(b)和文本文件(t)。在Windows下，打开文件的默认是二进制格式，如果要以文本方式打开，则必须在打方式中加上字符't'。

%fid = fopen('作业1.txt','wt+');%文件有两种格式：二进制文件(b)和文本文件(t)。在Windows 下，打开文件的默认是二进制格式，如果要以文本方式打开，则必须在打方式中加上字符't'。

fprintf(fid,'%c',' 班级姓名

学号作业序号主要内容');

fprintf(fid,'%c\n',' ');fprintf(fid,'%c',' 测绘12-1 王亚亚07122825 上机实验作业三精密单点定位');

fprintf(fid,'%c\n',' ');

for i=1:1%输出datee1矩阵即起始时间

for j=1:6

fprintf(fid,'%d\t',datee1(i,j));

if rem(i,1)==0 && rem(j,6)==0

fprintf(fid,'%c','%起始时间');

%fprintf(fid,'\r\n');

end

end

end

fprintf(fid,'%c\n',' ');

for i=1:1%输出datee矩阵即终止时间

for j=1:6

fprintf(fid,'%d\t',datee(i,j));

if rem(i,1)==0 && rem(j,6)==0

fprintf(fid,'%c','%终止时间');

%fprintf(fid,'\r\n');

end

end

end

fprintf(fid,'\n\nX: %12.3f Y: %12.3f

Z: %12.3f %s', me(1,1), me(2,1), me(3,1),'%观测站（wsrt0950）位置'); %输出X、Y、Z的平均值%12.3f这个是格式要求

fclose(fid);

print -deps easy3%保存图像为eps 格式

2、主要函数代码及其注释：recpo_ls函数代码：

function [pos, El, GDOP, basic_obs] =

recpo_ls(obs,sats,time,Eph)

% RECPO_LS 用伪距观测值最小二乘计算接收机位置。

% obs 该历元的卫星观测值

% sats obs对应的卫星的PRN号

% time 对应的卫星时刻

% Eph 星历信息

v_light = 299792458;

dtr = pi/180;

m = size(obs,1); % number of svs 卫星个数el = zeros(m,1);% identify ephemerides columns in Eph

for t = 1:m

col_Eph(t) = find_eph(Eph,sats(t),time); %在星历矩阵中，搜索卫星对应的列。

End

% 接收机位置和接收机钟差的初步求解

pos = zeros(4,1);

no_iterations = 6; %迭代6次

ps_corr = [];

sat_pos = [];

for iter = 1:no_iterations

A = [];

12

omc = []; % observed minus computed observation

for i = 1:m

k = col_Eph(i);

tx_RAW = time - obs(i)/v_light;

t0c = Eph(21,k);

dt = check_t(tx_RAW-t0c); % GPS时间超限或下溢的修复

tcorr = (Eph(2,k)*dt + Eph(20,k))*dt + Eph(19,k);

tx_GPS = tx_RAW-tcorr;

X = satpos(tx_GPS, Eph(:,k));% 用eph 给定的星历信息,计算t时刻的m卫星的X，Y，Z坐标

if iter == 1

traveltime = 0.072;

Rot_X = X;

trop = 0;

else

rho2 = (X(1)-pos(1))^2+(X(2)-pos(2))^2+(X(3)-pos(3))^2; %接收机与卫星间的距离

traveltime = sqrt(rho2)/v_light;

Rot_X =

e_r_corr(traveltime,X); % E_R_CORR返回卫星的ECEF坐标值，考虑信号传播时地球自转的影响

rho2 = (Rot_X(1)-

pos(1))^2+(Rot_X(2)-pos(2))^2+(Rot_X(3)-

pos(3))^2;

[az,el,dist] =

topocent(pos(1:3,:),Rot_X-

pos(1:3,:)); %topocent(X,dx)将向量dx转换成以X为坐标原点的站心坐标。

%两个参数都是3*1的向量

if iter == no_iterations, El(i) = el; end

trop =

tropo(sin(el*dtr),0.0,1013.0,293.0,50.0,...%T ROPO计算对流层改正

0.0,0.0,0.0);

end

% subtraction of pos(4) corrects for receiver clock offset

% pos(4) 接收机钟差改正

% v_light*tcorr is the satellite clock offset

% v_light*tcorr 卫星钟差改正

if iter == no_iterations

ps_corr = [ps_corr;

obs(i)+v_light*tcorr-trop];

sat_pos = [sat_pos; X'];

end

omc = [omc; obs(i)-norm(Rot_X-

pos(1:3),'fro')-pos(4)+v_light*tcorr-trop];

A = [A; (-(Rot_X(1)-pos(1)))/obs(i)...

(-(Rot_X(2)-pos(2)))/obs(i) ...

(-(Rot_X(3)-pos(3)))/obs(i) 1];

end% i

x = A\omc;

pos = pos+x;

if iter == no_iterations, GDOP =

sqrt(trace(inv(A'*A)));

%% two lines that solve an exercise on computing tdop

% invm = inv(A'*A);

% tdop = sqrt(invm(4,4))

end

end% iter

basic_obs = [sat_pos ps_corr];

tropo函数代码：

function ddr =

tropo(sinel,hsta,p,tkel,hum,hp,htkel,hhum)

%TROPO Calculation of tropospheric correction. 计算对流程改正

% The range correction ddr in m is to be subtracted from 在米级范围的ddr需要由伪距和载波相位做差得到

% pseudo-ranges and carrier phases % sinel sin of elevation angle of satellite 卫星高度角的正弦值

% hsta height of station in km 以千米为单位的测站高度（大地高？）

% p atmospheric pressure in mb at height hp 高度上升时候的大气压（单位mb）% tkel surface temperature in degrees Kelvin at height htkel 地表温度（单位凯氏度）以每

13

千米为测量单位下

% hum humidity in % at height hhum 湿度（单位%）在每千米为测量单位下

% hp height of pressure measurement in km 高大气压测量（单位千米）

% htkel height of temperature measurement in km 高温度测量(单位千米)

% hhum height of humidity measurement in km 高湿度测量（单位千米）

% tk 绝对温度

% hum 湿度

% tkel 温度测量单位km

% a 大气

% Goad, C.C. & Goodman, L. (1974) A Modified Tropospheric 改正后的对流层折射修正模型% A Matlab reimplementation of a C code from driver. 由驱动程序使得Matlab重新实现一个C 代码

a_e = 6378.137; % semi-major axis of earth ellipsoid 地球椭球的长半轴长度（单位千米）

b0 = 7.839257e-5; % 地球自转速率rad/s tlapse = -6.5; %大气温度下降率：高程每增加1km，温度下降6.5度

tkhum = tkel+tlapse*(hhum-htkel);%湿度影响下的绝对温度（指的是地面温度）=绝对温度+（-6.5°）*（湿大气层顶高-温度（单位km））

atkel = 7.5*(tkhum-273.15)/(237.3+tkhum-273.15); %atkel：（应用于下方e0的公式）大气影以及湿度影响下的数值（km）绝对零度，也就是-273.15℃(摄氏度)。绝对零度应比它低273.15℃所以精确的换算关系应该是T(K)热力学温标=t(℃)摄氏温标+273.15热力学温标和摄氏温标间的换算

e0 = 0.0611*hum*10^atkel; %地面水汽压：（mbar）参考：《三种对流层延迟改正模型精度评估》中第3页顶部

tksea = tkel-tlapse*htkel; %海面绝对温度=温度测量（km）-（-6.5°）*高温度测量（km）em = -978.77/(2.8704e6*tlapse*1.0e-5); %em：为对流层天顶方向的干延迟中需要计算的其中一项（幂次项）参考《EGNOS 对流层延迟改正模型及其精度分析》第一页公式（1）与公式（2）第一个g=9.806 65m/s2(重力加速度)化成厘米为单位：980.66、tlapse：大气温度下降率（每km下降转化为每cm下降，单位换算，乘以10的负5次方）2.8704e6 Rd=287.054 J/（kg. K）（K是凯氏温度）热量的单位是焦耳温度的单位是摄氏度不过,通过一个物体的"比热容"就可以转化了。比热容意思是,单位质量(1kg)的物体升高1摄氏度所要吸收的热量.

tkelh = tksea+tlapse*hhum;

e0sea = e0*(tksea/tkelh)^(4*em); %海面水汽压tkelp = tksea+tlapse*hp;

psea = p*(tksea/tkelp)^em;

if sinel < 0 %如果高度角正弦值小于0，修正为0

sinel = 0;

end;

tropo = 0;

done = 'FALSE';

refsea = 77.624e-6/tksea; %k1=77.604 K/mbar 用于计算对流层天顶方向的干延迟参考：《三种对流层延迟改正模型精度评估》

htop = 1.1385e-5/refsea;

refsea = refsea*psea; %计算对流层天顶方向的干延迟

ref = refsea*((htop-hsta)/htop)^4;

while 1

rtop = (a_e+htop)^2-(a_e+hsta)^2*(1-

sinel^2);%公式中Ri的前半部分

if rtop < 0, rtop = 0; end; % check to see if geometry is crazy

rtop = sqrt(rtop)-(a_e+hsta)*sinel;%公式中的Ri

a = -sinel/(htop-hsta);%公式中的a

b = -b0*(1-sinel^2)/(htop-hsta);%公式中的b

rn = zeros(8,1);

for i = 1:8

rn(i) = rtop^(i+1);

end;

alpha = [2*a, 2*a^2+4*b/3, a*(a^2+3*b),...

a^4/5+2.4*a^2*b+1.2*b^2,

2*a*b*(a^2+3*b)/3,...

b^2*(6*a^2+4*b)*1.428571e-1, 0, 0];

14

if b^2 > 1.0e-35, alpha(7) = a*b^3/2; alpha(8) = b^4/9; end;

dr = rtop;

dr = dr+alpha*rn;

tropo = tropo+dr*ref*1000;;%干延迟与湿延迟的和

if done == 'TRUE ', ddr = tropo; break; end;

done = 'TRUE ';

refsea = (371900.0e-6/tksea-12.92e-6)/tksea;

htop = 1.1385e-

5*(1255/tksea+0.05)/refsea;;%湿对流层高度ref = refsea*e0sea*((htop-hsta)/htop)^4;;%对流层湿分量折射系数

end;

end

四、实验总结：

通过编程实现GPS单点定位，我们更加深入的了解到星历数据的各个意义，而且通过此次试验加深了我们对编程的认识和对MATLAB的认识。我们认识到MATLAB 编程可以减轻我们大量的工作负担。我们也深刻的了解到学校为把我们培养成出色的工程师对我们的严要求。通过此次试验也让我们深刻体会到测绘专业是一门综合型专业，我们要具有很强的编程能力，要注重对编程的培养，这样才能在以后的工作中更好的处理工作遇到的问题。为以后从事测绘工作奠定了一定的基础。

15