卫星导航GPS典型例题编程报告

GPS Matlab编程报告

一、通过星历计算卫星坐标（修正） 1、 算法流程

（1） 计算GPS卫星运行的平均速度n

n?n0??n n0?（2） 计算归化时间?t

?a3,a??a?

2s?t?t?toe

（3） 计算计算观测历元t的平近点角M

M?M0?n?t

（4） 计算偏近点角E

利用不动点迭代法求解此方程：E?M?esinE 得出E 不动点迭代法的迭代公式为：xk?1???xk? k?0,1,2?

（5） 计算卫星的地心矢径r0

r0?a?1?ecosE?

（6） 计算真近点角f

f?2*arctan(（7） 计算升交点角距?0

1?e?E?tan??) 1?e?2??0???f

（8） 计算摄动改正项：?u,?r,?i

?u?Cussin2?0?Cuccos2?0?r?Crssin2?0?Crccos2?0 ?i?Cissin2?0?Ciccos2?0（9） 计算经过摄动改正的升交点角距?，卫星矢径r，和轨道面倾角i

???0??ur?r0??r??ti?i0???i（10） 计算观测历元t的升交点经度?k

?k??0?????ie??t??ietoe

（11） 计算卫星在轨道平面坐标系中的位置?x0,y0,z0?

?x0??cos??????y?rsin??0??? ?z?? 0??0???（12） 计算卫星在地固坐标系下的坐标?x,y,z?

?x??x0??x0?cos?k?y0?cosi?sin?k???????y?R(??)R(?i)y?x?sin??y?cosi?cos?zkx00k0k?????? ?z??z???y0?sini?z0cosi???0???

2、 Matlab源程序

%本程序用于通过星历计算卫星坐标（修正） %部分星历数据未用到 clear;

%---------t时刻，卫星星历数据---------

t=1.728128741751984e+005; %当前时间（接收到卫星信号时间） Omega0= -3.1122016819; %按参考时间计算的升交点经度 I0=0.3*pi; %非GEO卫星 %I0= 0.0879969068; %GEO卫星 SqrtA= 6493.27499625; %长半轴平方根 Ecc= 0.01995786; %偏心率 Small_Omega= -1.7409228484; %近地点角距 Mu0= 0.0005884201; %参考时间的平近点角

Delta_n= -2.4123e-009; %卫星平均运动速度与计算值之差 I_Dot= -3.4333e-010; %轨道倾角变化率 Omega_dot= 3.47676e-009; %升交点经度变化率

C_uc= 4.692e-006; %纬度幅角的余弦调和改正项的振幅 C_us= -1.10036e-005; %纬度幅角的正弦调和改正项的振幅 C_ic= 9.01e-008; %轨道倾角的余弦调和改正项的振幅 C_is= 2.01e-008; %轨道倾角的正弦调和改正项的振幅 C_rc= 325.73368; %轨道半径的余弦调和改正项的振幅 C_rs= 146.91074; %轨道半径的正弦调和改正项的振幅

Toe= 172800; %星历参考时间 IODC= 10; %钟差数据龄期 URAI= 0; %用户距离精度标志 IODE= 10; %星历数据龄期

Toc = 172800; %本时段钟差参数参考时间 a0 = -2.72893e-005; %卫星钟差改正0阶多项式系数 a1 = -6.7531e-014; %卫星钟差改正1阶多项式系数 a2 = 5.787e-018; %卫星钟差改正2阶多项式系数

%---------定义常量--------- c=2.99792458e8; %光速

omegae=7.2921151467e-5; %地球自转角速度 mu=3.986004418e14; %地球引力常数GM

%---------1、计算GPS卫星运行的平均速度n--------- a=SqrtA*SqrtA;

n=sqrt(mu/(a^3))+Delta_n;

%---------2、计算归化时间Delta_t--------- Delta_t=t-Toe;

%---------3、计算观测历元t的平近点角M--------- M=Mu0+n*Delta_t;

%---------4、计算偏近点角E--------- eps=1e-20; E=M; tol=1;

while (tol>eps) %不动点迭代法 E0=E;

E=M+(Ecc)*sin(E0); tol=abs(E-E0); end

%---------5、计算卫星的地心矢径r0--------- r0=a*(1-Ecc*cos(E));

%---------6、计算真近点角f--------- %f=2*atan(sqrt((1+Ecc)/(1-Ecc))*tan(E/2)); f=atan((sqrt(1-Ecc^2)*sin(E))/(cos(E)-Ecc));

%---------7、计算升交点角距Phi0--------- Phi0=Small_Omega+f;

%---------8、计算摄动改正项：Delta_u,Delta_r,Delta_i--------- Delta_u=C_us*sin(2*Phi0)+C_uc*cos(2*Phi0); Delta_r=C_rs*sin(2*Phi0)+C_rc*cos(2*Phi0); Delta_i=C_is*sin(2*Phi0)+C_ic*cos(2*Phi0);

%---------9、计算经过摄动改正的升交点角距Phi，卫星矢径r，和轨道面倾角I--------- Phi=Phi0+Delta_u; r=r0+Delta_r;

I=I0+Delta_i+I_Dot*Delta_t;

%---------10、计算观测历元t的升交点经度Omegak--------- Omegak=Omega0+(Omega_dot-omegae)*Delta_t-omegae*Toe;

%---------11、计算卫星在轨道平面坐标系中的位置--------- x0=r*cos(Phi); y0=r*sin(Phi); z0=0;

%---------12、计算卫星在地固坐标系下的坐标--------- x=x0*cos(Omegak)-y0*cos(I)*sin(Omegak); y=x0*sin(Omegak)+y0*cos(I)*cos(Omegak); z=y0*sin(I)+z0*cos(I);

%---------输出卫星坐标---------

fprintf ('（修正后）卫星在地固坐标系中的坐标：\\n X=%.10f Y=%.10f Z=%.10f\\n',x,y,z);

%老师给的结果：卫星位置：X=7073881.4181256806 Y=23901970.8378255780 Z=-32955601.5025106560 X=7073881.4181256806;Y=23901970.837825578;Z=-32955601.5025106560;%修正后的，非GEO fprintf ('与老师的结果的偏差为：\\n D_X=%.10f D_Y=%.10f D_Z=%.10f\\n',x-X,y-Y,z-Z);

3、 程序运行结果

（修正后）卫星在地固坐标系中的坐标：

X=7073887.5144389411 Y=23901969.0496121940 Z=-32955601.4908931700 与老师的结果的偏差为：

D_X=6.0963132605 D_Y=-1.7882133834 D_Z=0.0116174854

二、已知4颗卫星坐标及测得的伪距（已改正）求接收机位置 1、 原理及算法

不考虑电离层延迟和对流层延迟的因素时，由4颗卫星在地心坐标系中的坐标及对应的伪距即可计算出接收机的位置，按如下四元二次方程组求解四个未知数：

Ri?R?c??t0??xi?x???yi?y???zi?z?222?c??t0 i?1,2,3,4

其中，Ri为已改正的伪距，R为卫星与接收机之间实际的距离，(xi,yi,zi)分别为4颗卫星的坐标，?x,y,z?为接收机在地心坐标系中的坐标，?t0为接收机时钟与卫星时钟的偏差值。

采用Newton迭代法求解此四元二次非线性方程组，其算法为： 将上述方程组写成如下形式四元方程组：

?f1?x,y,z,?t0??0??f2?x,y,z,?t0??0 ??f3?x,y,z,?t0??0?fx,y,z,?t?00??4?记向量X??x,y,z,?t0?，F??f1,f2,f3,f4?，则上式化为：F?X??0 若给出此方程组的一个初值X?k?TT?x,y,z,?t0??k??k??k??k??T，将函数F?X?的分量fi?X?在

X?k?处Taylor展开，并取线性部分，则可表示为：

F?X??FX?k??F?X?k?X?X?k?

再由F?X??0，可认为：

其中Jacobi矩阵为：

???????F?X?k?X?X?k???FX?k?

?????（1）

??f1??x???f2??xF??X?????f3???x??f4???x则（1）方程组的解为：

?f1?y?f2?y?f3?y?f4?y?f1?z?f2?z?f3?z?f4?z?f1???t0???f2???t0?? ?f3????t0??f4????t0?

本文来源：https://www.bwwdw.com/article/xip7.html