---

前言

为了研究GNSS卫星的广播星历的精度，需要获得GNSS卫星的坐标与实际位置的差值。使用MATLAB编程，实现提取广播星历参数和精密星历精确坐标的过程，并基于利用广播星历计算GNSS卫星位置的原理计算得到GNSS卫星的瞬时坐标。将计算得到的GNSS卫星瞬时坐标与编程提取到的同一时刻的精密星历中的精确坐标对比进行精度分析，进而得到GNSS卫星广播星历的精度。

提示：以下是本篇文章正文内容，下面案例可供参考

一、GNSS卫星坐标计算的程序实现

（1）利用stuct（）函数建立空结构体，用以存储读取导航电文获取的参数，根据卫星计算原理的不同，结构体存储变量的数目也不同。
（2）从头文件结束标志‘END OF HEADER’开始，遍历存储导航电文，将每个卫星的参数赋给相应的结构体，同时使用str2num（）函数将字符变为数值，方便后续计算。
（3）选取要计算的卫星，调用对应结构体里面的参数，作为输入参数赋给即将调用的函数 。函数，从而获得返回值，即卫星的坐标。可调用的函数有GPScucalation（），bdscaculation（），glonasscaculation（）。
(4)将计算的卫星坐标统一放在一个结构体中，便于后续的精度分析。

1. 基于 Matlab 读取标准 RINEX 格式星历数据

计算GNSS卫星的瞬时位置需要用到卫星的广播星历，而广播星历多为RINEX格式。想要获得所需要的参数，还需要对RINEX 格式进一步处理。而使用matlab对RINEX数据进行处理的思路，是根据RINEX格式本身的特点，从特殊字符串‘END OF HEADER’头文件的结束标志开始，通过简单的循环语句，将星历数据存储进预先定义好的结构体中。具体代码实现如下：
（1）通过stuct（）函数创建存储星历数据的空结构体

navdata=struct('prn',nan,'year',nan,'month',nan,'day',nan,'hour',nan,'min',nan,'sec',nan,'a0',nan,'a1',nan,'a2',nan,'aode',NaN,'crs',NaN,'dn',NaN,.'m0',NaN,'cuc',NaN,'e',NaN,'cus',NaN,'a',NaN,'toe'NaN,'cic',NaN,'om0',NaN,'cis',NaN,'i0',NaN,'crc',NaN,'w',NaN,'dom',NaN,'di',NaN,'cflg12',NaN,'weekno',NaN,'pflg12',NaN,'svaccuracy',NaN,'svhealth',NaN,'tgd',NaN,'aodc',NaN,'transmission', NaN,'s1',NaN,'s2',NaN,'s3',NaN);

所创建结构体的字段名与星历提供的参数有关，GPS卫星，GALILEO卫星与北斗卫星都可以采用这段函数来存放他们自己的星历数据，字符‘nan’是不为该字符指定任何一种格式。而glonass因为星历所提供的参数不同，创建结构体所用的字段名自然也不同，其创建结构体的代码如下：

navdata=struct('prn',nan,'year',nan,'month',nan,'day',nan,'hour',nan,'min',nan,'sec',nan,'sv1',nan,'sv2',nan,'utc',nan,'x',NaN,'u',NaN,'ax',NaN,'bn',NaN,'y',NaN,'v',NaN,'ay',NaN,'icd',NaN,'z',NaN,'w',NaN,'az',NaN,'age',NaN);

（2）打开导航电文后，通过简单的while循环，利用findstr（）函数逐行查找特殊字符‘END OF HEADER’，获得读取数据的起始位置。因为GPS与GLONASS卫星的广播星历所采用的为rinex2.02的格式，所以可以在获取函数位置的同时，可以同时获取星历的行数，来计算星历中卫星的个数。GALILEO卫星与北斗卫星采用的为rinex3.02格式，可直接通过findstr（）函数来获取起始位置，部分代码如下。

fid = fopen('rinex');
while 1 
headlines = headlines+1;
line = fgetl(fid);
answer = findstr(line,'END OF HEADER');
if ~isempty(answer)
break;
end
end
while 1
noeph = noeph+1;
line = fgetl(fid);
if line == -1
break
end
end
noeph = noeph/8;

在计算glonass卫星个数时，因为星历文件只提供四行参数，所以应使用noeph = noeph/4语句。
（3）通过for循环使空结构体数组中能够存储获取GNSS卫星的星历数据。最后能通过简单的循环语句来读取调用各卫星星历数据。对于GPS与glonass卫星所采用的的rinex2.0格式星历文件，可以通过计算卫星数来建立一个简单的for循环，从而逐行获取数据，每8行或者4作为一组数据。而对于北斗与GALILEO卫星所采用的rinex3.02格式可通过findstr（）函数匹配字符串寻找所需要的卫星参数，同样每8行可以作为一组数据。星历每列对应的数据作为提取的依据。部分代码如下：

for l = 1:head_lines
line = fgetl(fid);
end
for i = 1:noeph
line = fgetl(fid); 
navdata(i).prn = str2num(line(1:2));
y = str2num(line(3:6));
if y >79
navdata(i).year = y+1900;
else
navdata(i).year = y+2000;
end
navdata(i).m = str2num(line(7:9));
navdata(i).d = str2num(line(10:12));
navdata(i).h = str2num(line(13:15));
navdata(i).min = str2num(line(16:18));
navdata(i).sec = str2num(line(19:22));
end
status = fclose(fid);

同样注意，由于rinex2.0格式与rinex3.02格式不同，所以在使用line（）函数提取相应的字符时应注意是否对应参数所在的列数。同时需要使用str2num（）函数将字符型转化为数值型，提取的参数才能正常使用。另外还需要注意，在提取GLONASS函数时需要注意，因为参数的单位为千米，所以要将其化为米，方便后续计算。
(4) 要想分析广播星历的精度，还要获得精密星历对应时刻的准确坐标值来作比较。精密星历的读取与广播星历的读取类似，都是通过简单的循环结构依次给空结构体赋值，但要注意的是精密星历读取不用寻找开始位置，直接跳过32行开始赋值。另外精密星历的单位为千米，需要转换为米，其实现的部分代码如下：

for i=1:32 line=fgetl(fid);
end
while 1 line=fgetl(fid);result=findstr(line,'EOF');if(isempty(result))SP3_data(t).year=str2num(line(4:7));for j=1: nasdfgmumline=fgetl(fid);if ~isempty(findstr(line,'PG'))SP3_data(t).sate(j)=0;endSP3_data(t).PRN(j)=str2num(line(3:4));SP3_data(t).x(j)=str2num(line(5:18)); end
end
fclose(fid);

nasdfgmum为该星历所含卫星的个数，每个星历含有的卫星数目不同。另外根据需要将每种卫星赋予标记值，GPS卫星的标记为0, GLONASS卫星的标记为1,Galileo卫星的标记为2， BDS卫星的标记为3，其余卫星的标记为4.
想要获得卫星的精密星历的坐标，可以遍历结构体，通过判定给每种卫星赋予的特定的标记，可以获得特定卫星在24h内的全部坐标。

for i=1:tfor k=1:nasdfgmumif SP3_data(i).minute==0if SP3_data(i).sate(k)==0&&SP3_data(i).PRN(k)==01temp(m,1)= SP3_data(i).x(k)*1000;temp(m,2)=SP3_data(i).y(k)*1000;temp(m,3)=SP3_data(i).z(k)*1000;m=m+1;endendend
end

2基于Matlab的GPS/GALILEO卫星位置计算

根据基于广播星历计算卫星位置的原理来设计matlab程序，同时在使用matlab计算GPS卫星的位置时，需要注意一些问题：
（1）计算真近点角时，要使用四象限反正切函数atan2（）来计算。
（2）偏近点角的计算需要进行迭代，使用while…end函数，当满足 ∣Ek+1-Ek∣< 10-12时，停止循环。
GALILEO卫星位置的计算原理与GPS相同，所以此函数同样可以应用计算GALILEO卫星的位置。虽然galieo卫星所采用的坐标系与GPS略有不同，但因为影响很小，所以可以忽略掉此项误差，GPS卫星与GALILEO卫星因为采用时间系统不同而造成的坐标间的误差可以忽略。
GPS/GALILEO卫星位置计算函数gpsposition部分代码如下：

n0=sqrt(GM)/(a^3); 
n=n0+dn;
UT=hour+(min/60)+(sec/3600); 
julianday=fix(362.25*year)+fix(30.6001*(month+1))+day+(UT/24)+1720981.5; 
t=fix((julianday -2444244.5)/7); 
t=( julianday -2444244.5-t*7)*24*3600; 
t=t-(a0+a1*(t-toe)+a2*(t-toe)*(t-toe)); 
tk=t-toe;
if tk>302400tk=tk-604800;
elseif tk<-302400tk=tk+604800;
endmk=m0+n*tk; ek=mk;ek0=inf;while abs(ek0-ek)>=1e-12ek0=ek;ek=mk+e*sin(ek0);endsinfk=(sqrt(1-e*e)*sin(ek))/(1-e*cos(ek));  cosfk=(cos(ek)-e)/(1-e*cos(ek));fk=atan2(sinvk,cosvk);faik=vk+w;deltu=cuc*cos(2*faik)+cus*sin(2*faik); deltr=crc*cos(2*faik)+crs*sin(2*faik);delti=cic*cos(2*faik)+cis*sin(2*faik);uk=faik+deltu; rk=a*a*(1-e*cos(ek))+deltr;ik=i0+delti+di*tk;
xk=rk*cos(uk); 
yk=rk*sin(uk);
zk=0;
lk=om0+(dom-we)*tk-we*toe; 
X=xk*cos(lk)-yk*cos(ik)*sin(lk) 
Y=xk*sin(lk)+yk*cos(ik)*cos(lk)
Z=yk*sin(ik

要调用此函数时，可以使用代码：

position=gpsposition( navdata(i).y,navdata(i).m,navdata(i).d,navdata(i).h,navdata(i).min,navdata(i).sec,navdata(i).a0,navdata(i).a1,navdata(i).a2,navdata(i).crs,navdata(i).dn,navdata(i).m0,navdata(i).cuc,navdata(i).e,navdata(i).cus,navdata(i).a,navdata(i).toe,navdata(i).cic,navdata(i).om0,navdata(i).cis,navdata(i).i0,navdata(i).crc,navdata(i).w,navdata(i).dom,navdata(i).di);

3 基于Matlab的北斗卫星位置计算

由于3种不同的卫星系统组成了北斗卫星导航系统，其中IGSO卫星与MEO卫星的matlab计算程序可以用计算GPS卫星坐标的原理计算，而GEO卫星位置的计算可以根据第二章的原理利用以下程序计算。使用matlab计算北斗卫星时应注意：
（1）在计算归化时间时，需要将民用日转化为儒略日，而因为北斗卫星与GPS卫星的时间系统并不统一，需要将归化时间减去14s来将北斗卫星的时间系统纳入到GPS时间系统中。
（2）在计算Rx旋转矩阵时，因为-5°不是弧度，需要使用cosd（）/sind（）函数计算正余弦的值。
（3）因为要计算三种卫星,所以在计算卫星的位置时，要使用if…elseif…elsef…else…end函数
北斗卫星GEO卫星部分计算代码如下:

lk=om0+(dom)*tk-we*toe;
%计算BDS卫星在坐标系中的坐标
Xg=xk*cos(lk)-yk*cos(ik)*sin(lk);
Yg=xk*sin(lk)+yk*cos(ik)*cos(lk);
Zg=yk*sin(ik);
temp=[Xg;Yg;Zg];
Rx=[1,0,0;0,cosd(-5),sind(-5);0,-sind(-5),cosd(-5)];
Rz=[cos(we*tk),sin(we*tk),0;-sin(we*tk),cos(we*tk),0;0,0,1];
temp=Rz*Rx*temp;
X=temp(1);
Y=temp(2);
Z=temp(3);
result=[X,Y,Z]; 

调用此函数的方法与GPS卫星计算坐标函数的调用方式相同

4 基于Matlab的GLONASS卫星位置计算

GLONASS卫星的计算原理与其他卫星都不同，大多数卫星的瞬时坐标是根据参考时刻的开普勒轨道根数与摄动改正数按照一定的原理计算的，而由于GLONASS卫星星历只提供卫星参考时刻的三维坐标，速度向量，以及日月摄动加速度。所以要计算GLONASS卫星任意时刻的位置，需要采用轨道数值积分的方法。一般采用经典四阶龙贝格-库塔方法来计算GLONASS卫星位置。
计算GLONASS卫星坐标需要注意的地方是，因为GPS时与GLONASS时的时间系统起点不同，而GLONSS时跟着UTC时间一起变化所以存在跳秒现象[7]。所以在计算积分长度时需要减去17s来将GLONASS时并入到GPS时中。同时因为坐标系统不同带来的误差很小，所以在具体计算时不计入误差。考虑积分出来的精度问题，积分步长选择30s，积分长度选择30min[5]。
计算glonass卫星坐标的部分代码为：

function result = gloposition( x,y,z,u,v,w,xls,yls,zls )
length =30*60-17;
h=30;
for i=h:h:lengtha1=f(x,y,z,u,v,w,xls,yls,zls);u1=u;v1=v;w1=w;    
a2=f(x+u1*h/2,y+v1*h/2,z+w1*h/2,u+a1(1)*h/2,v+a1(2)*h/2,w+a1(3)*h/2,xls,yls,zls);u2=u+a1(1)*h/2;v2=v+a1(2)*h/2;w2=w+a1(3)*h/2;   a3=f(x+u2*h/2,y+v2*h/2,z+w2*h/2,u+a2(1)*h/2,v+a2(2)*h/2,w+a2(3)*h/2,xls,yls,zls);u3=u+a2(1)*h/2;v3=v+a2(2)*h/2;w3=w+a2(3)*h/2;    
a4=f(x+u3*h,y+v3*h,z+w3*h,u+a3(1)*h,v+a3(2)*h,w+a3(3)*h,xls,yls,zls);u4=u+a3(1)*h;v4=v+a3(2)*h;w4=w+a3(3)*h;u=u+h*(a1(1)+2*a2(1)+2*a3(1)+a4(1))/6;v=v+h*(a1(2)+2*a2(2)+2*a3(2)+a4(2))/6;w=w+h*(a1(3)+2*a2(3)+2*a3(3)+a4(3))/6;x=x+h*(u1+2*u2+2*u3+u4)/6;y=y+h*(v1+2*v2+2*v3+v4)/6;z=z+h*(w1+2*w2+2*w3+w4)/6; 
end
result=[x,y,z]

其中，所调用的f（）函数的部分代码如下：

format long e;
r = sqrt(xi*xi+yi*yi+zi*zi);
ax=-GM*xi/(r^3)-1.5*C02*GM*ae*ae*xi*(1-(5*(zi^2))/(r^2))/(r^5)+we*we*xi+2*we*vy+xls;
ay=-GM*yi/(r^3)-1.5*C02*GM*ae*ae*yi*(1-(5*(zi^2))/(r^2))/(r^5)+we*we*yi-2*we*vx+yls;
az = -GM*zi/(r^3)-1.5*C02*GM*ae*ae*zi*(3-(5*(zi^2))/(r^2))/(r^5)+zls;
result = [ax,ay,az];

调用此函数时，可使用以下代码：

position=gloposition( navdata(i).x,navdata(i).y,navdata(i).z,navdata(i).u,navdata(i).v,navdata(i).w,navdata(i).ax,navdata(i).ay,navdata(i).az);

计算GLONASS卫星位置时，也可以使用MATLAB的内置函数ODE45来计算，与通过经典四阶龙贝格-库塔计算出来的结果相差不大。

总结

全部代码在我上传的资源里面，可以自行下载，代码运行的一些细节在这篇文章里也有，这篇文章仅做交流参考用。