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

GPSMatlab编程报告一、通过星历计算卫星坐标（修正）1、算法流程（1）计算GPS卫星运行的平均速度n0nnn203,snaaa（2）计算归化时间toettt（3）计算计算观测历元t的平近点角M0MMnt（4）计算偏近点角E利用不动点迭代法求解此方程：sinEMeE得出E不动点迭代法的迭代公式为：10,1,2kkxxk（5）计算卫星的地心矢径0r01cosraeE（6）计算真近点角f12*arctan(tan)12eEfe（7）计算升交点角距00f（8）计算摄动改正项：,,uri000000sin2cos2sin2cos2sin2cos2uusucrrsrciisicCCCCCC（9）计算经过摄动改正的升交点角距，卫星矢径r，和轨道面倾角i000urrriiit（10）计算观测历元t的升交点经度k0kieieoett（11）计算卫星在轨道平面坐标系中的位置000,,xyz000cossin0xyrz（12）计算卫星在地固坐标系下的坐标,,xyz000000000coscossin()()sincoscossincoskkzkxkkxxxyiyRRiyxyizzyizi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(toleps)%不动点迭代法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('（修正后）卫星在地固坐标系中的坐标：\nX=%.10fY=%.10fZ=%.10f\n',x,y,z);%老师给的结果：卫星位置：X=7073881.4181256806Y=23901970.8378255780Z=-32955601.5025106560X=7073881.4181256806;Y=23901970.837825578;Z=-32955601.5025106560;%修正后的，非GEOfprintf('与老师的结果的偏差为：\nD_X=%.10fD_Y=%.10fD_Z=%.10f\n',x-X,y-Y,z-Z);3、程序运行结果（修正后）卫星在地固坐标系中的坐标：X=7073887.5144389411Y=23901969.0496121940Z=-32955601.4908931700与老师的结果的偏差为：D_X=6.0963132605D_Y=-1.7882133834D_Z=0.0116174854二、已知4颗卫星坐标及测得的伪距（已改正）求接收机位置1、原理及算法不考虑电离层延迟和对流层延迟的因素时，由4颗卫星在地心坐标系中的坐标及对应的伪距即可计算出接收机的位置，按如下四元二次方程组求解四个未知数：222001,2,3,4iiiiRRctxxyyzzcti其中，iR为已改正的伪距，R为卫星与接收机之间实际的距离，(,,)iiixyz分别为4颗卫星的坐标，,,xyz为接收机在地心坐标系中的坐标，0t为接收机时钟与卫星时钟的偏差值。采用Newton迭代法求解此四元二次非线性方程组，其算法为：将上述方程组写成如下形式四元方程组：10203040,,,0,,,0,,,0,,,0fxyztfxyztfxyztfxyzt记向量0,,,TxyztX，1234,,,TffffF，则上式化为：FX0若给出此方程组的一个初值0,,,TkkkkkxyztX，将函数FX的分量ifX在kX处Taylor展开，并取线性部分，则可表示为：kkkFXFXFXXX再由FX0，可认为：kkkFXXXFX（1）其中Jacobi矩阵为：11110222203333044440ffffxyztffffxyztffffxyztffffxyztFX则（1）方程组的解为：11=0,1,2,kkkkkXXFXFX（2）利用（2）式进行迭代运算即可解出满足一定精度要求的近似解。若考虑地球自转影响，需要确定地球的自转角速度e以及发射信号瞬时到接收信号瞬时经历的时间t，从而计算出卫星信号传输时间内地球自转过的角度et，然后利用旋转矩阵I将上述计算得到的坐标值,,xyz做一旋转变换，得出考虑地球自转影响时接收机位置,,rrrxyz，其旋转操作如下：cossin0sincos0001rrrxxxyyyzzzI其中卫星信号传输总时间t可由最大伪距除以光速得到：max/itRc，（c为光速）。2、Matlab源程序functionReceiver_Position%本程序求解已知4颗卫星坐标及测得的伪距时接收机的位置坐标[XYZ]，需要求解四元二次非线性方程组%利用Newton法求解此非线性方程组clear;tic;c=2.99792458e8;%光速omegae=7.2921151467e-5;%地球自转角速度x0=[0;0;0;0];%取初始值tol=1.0e-6;%计算精度x1=x0-inv(df(x0))*fc(x0);%计算第一个迭代值n=1;while((norm(x1-x0)=tol)&&n10^8)%牛顿迭代计算x0=x1;x1=x0-inv(df(x0))*fc(x0);%迭代公式n=n+1;endT=toc;fprintf('接收机坐标为：X=%.10fY=%.10fZ=%.10f\n',x1(1),x1(2),x1(3));%老师的计算结果：x=1725670.7674292959y=-2116958.3717429629z=3129817.7967605507x=1725670.7674292959;y=-2116958.3717429629;z=3129817.7967605507;fprintf('与老师的结果的偏差为：D_X=%.10fD_Y=%.10fD_Z=%.10f\n',x1(1)-x,x1(2)-y,x1(3)-z);fprintf('迭代次数为：n=%i\n',n);fprintf('计算时间为：T=%fs\n',T);fprintf('计算精度为：tol=%f\n\n',tol);%考虑地球自转因素R=[2.4310764064e+007,2.2914600784e+007,2.0628809405e+007,2.3422377972e+007];%伪距D_t=max(R)/c;%取最大的伪距除以光速得到卫星信号传输时间D_phi=omegae*D_t;%发射信号瞬时到接收信号瞬时，地球偏转的角度I=[cos(D_phi),sin(D_phi),0;-sin(D_phi),cos(D_phi),0;0,0,1];%旋转矩