作业报告--卫星坐标计算

-1-GPS卫星导航定位技术与方法作业报告之卫星坐标计算-2-1作业任务------------------------------------------------------------------------------------32作业思想---------------------------------------------------------------------------------------33作业条件及数据--------------------------------------------------------------------34作业过程---------------------------------------------------------------------------45源程序-----------------------------------------------------------------------------56计算结果---------------------------------------------------------------------------97心得体会与建议-----------------------------------------------------------------------------9-3-1作业任务根据教材所给广播星历参数，编程计算UTC2004年1月30日8点0分00秒—20分00秒，每隔1分钟的卫星号为7的卫星坐标。2作业思想根据教材提供广播星历参数，按照卫星坐标计算步骤一步一步求解表示卫星位置的量，最后求出卫星在地固坐标系中的空间直角坐标。根据历元间隔，计算出共21个历元的卫星坐标。此次作业的已知条件及计算步骤均为教材提供且思路清晰，简明易懂，只要在编程过程中注意个已知量、中间量和待求量的表示及数据对应正确，涉及角度的计算时小心避免错误即可。3作业条件及数据由GPS卫星导航文件（表3.1）提取出7号卫星的轨道参数及其摄动改正项表3.1提取及整理后的广播星历参数：-4-星历参数参数值参数含义Δn4.85305928777E-09给定的星历参数计算值的平均速度no之改正数(弧度)√a5.15373039818E+03卫星轨道长半轴的平方根t0e4.60800000000E+05星历的参考时刻，从星期日零时开始量度，以秒计M0-2.14992476872E+00参考时刻t0e的平近点角（弧度）e1.25958691351E-02轨道偏心率ω-1.85137234909E+00近地点角距（弧度）Cuc2.75485217571E-06升交点角距的余弦调和改正项的振幅(弧度)Cus9.03755426407E-06升交点角距的正弦调和改正项的振幅(弧度)Crc1.91875000000E+02轨道半径的余弦调和改正项的振幅（弧度）Crs5.37812500000E+01轨道半径的正弦调和改正项的振幅（弧度）Cic2.06753611565E-07轨道倾角的余弦调和改正项的振幅（弧度）Cis-7.26431608200E-08轨道倾角的正弦调和改正项的振幅（弧度）IDOT-1.31434044204E-10轨道倾角的变率i09.37690094716E-01参考时刻的轨道倾角（弧度）Ω-7.91890109042E-09升交点赤经的变率（弧度/s）Ω01.93198315097E+00GPS周开始时刻的升交点经度（弧度）4作业过程-5-（10）升交点经度：λt=Ω0+（Ω-ωe）（t-t0e）-ωet0e式中：ωe为地球自转角速度，ωe=7.2921151467×10-5rad/s；Ω为升交点赤经变化率，Ω0为GPS周开始时刻的升交点经度，Ω和Ω0均由星历文件给出；t0为参考时刻，一般取t0=0，为一周的开始（星期日的子夜）；t0e为广播星历的参考历元时刻，从t0开始起算；（11）卫星在地固坐标系中的空间直角坐标为：()()0kkkkkkkXxYRRiyZ式中：()kR，()kRi为旋转矩阵，将其代入并展开后得coscossinsincoscossinkkkkkkkkkkkkkkkXxyiYxyiZyi5源程序#includeiomanip.h#includestdlib.h#includemath.h#includefstream.h#includeiostream.hconstintk=20;//历元数doubleX[k],Y[k],Z[k];//计算所得卫星坐标constdoubleGM=3986004.418E8;//地球引力常数constdoublewe=7.2921151467E-5;//地球自转角速度doubledn,a05,t0e,M0,e,w,Cuc,Cus,Crc,Crs,Cic,Cis,IDOT,i0,WMG,WMG0;//定义广播星历参数//广播星历参数导入voidInput(){ifstreamf1(广播星历参数.txt);if(!f1){cerr广播星历参数.txtfilenotopen!endl;-6-exit(1);}f1dn;f1a05;f1t0e;f1M0;f1e;f1w;f1Cuc;f1Cus;f1Crc;f1Crs;f1Cic;f1Cis;f1IDOT;f1i0;f1WMG;f1WMG0;f1.close();}//卫星坐标计算voidCompute(doublet,intk){doublen0,n,tk,Mk,E0,Ek,Vk,fk,du,dr,di,uk,rk,ik,xk,yk,Lt;//定义卫星坐标计算过程中间量//（）计算平均角速度n：n0=sqrt(GM/pow(a05,6));n=n0+dn;//（）计算规化时刻tk：tk=t-t0e;//（）计算平近点角Mk：Mk=M0+n*tk;//（）迭代计算偏近点角Ek：E0=Mk;do-7-{Ek=Mk+e*sin(E0);E0=Ek;}while(Mk+e*sin(E0)-E0=1e-8);//（）计算真近点角Vk：doublev1=sqrt(1-e*e)*sin(Ek);doublev2=cos(Ek)-e;Vk=atan2(v1,v2);//（）计算升交角距fk：fk=Vk+w;//（）摄动改正：du=Cuc*cos(2*fk)+Cus*sin(2*fk);//升交角距改正dudr=Crc*cos(2*fk)+Crs*sin(2*fk);//轨道向径改正drdi=Cic*cos(2*fk)+Cis*sin(2*fk);//轨道倾角改正diuk=fk+du;//改正后的升交角距ukrk=a05*a05*(1-e*cos(Ek))+dr;//改正后的轨道向径rkik=i0+di+IDOT*tk;//改正后的轨道倾角ik//（）计算卫星在升交点轨道直角坐标系的坐标：xk=rk*cos(uk);yk=rk*sin(uk);//（）计算升交点精度：Lt=WMG0+(WMG-we)*(t-t0e)-we*t0e;//（）计算卫星空间直角坐标：X[k]=xk*cos(Lt)-yk*cos(ik)*sin(Lt);Y[k]=xk*sin(Lt)+yk*cos(ik)*cos(Lt);Z[k]=yk*sin(ik);}voidmain()-8-{Input();for(inti=0;ik;i++){doublet=t0e+i*60;Compute(t,i);}cout历元/s卫星号X/mY/mZ/m卫地距离endl;for(i=0;ik;i++){cout.precision(11);coutt0e+i*607X[i]Y[i]Z[i]sqrt(X[i]*X[i]+Y[i]*Y[i]+Z[i]*Z[i])endl;}//输出结果到“卫星坐标计算结果.txt”ofstreamoutfile;outfile.open(卫星坐标计算结果.txt);if(outfile.is_open()){outfile历元/s卫星号X/mY/mZ/m卫地距离endl;for(inti=0;ik;i++){outfile.precision(11);outfilet0e+i*607X[i]Y[i]Z[i]sqrt(X[i]*X[i]+Y[i]*Y[i]+Z[i]*Z[i])endl;}}outfile.close();}-9-6计算结果表6.1卫星坐标计算结果历元/s卫星号X/mY/mZ/m4608007-13413856.78016092624.86716628017.8594608607-13535898.78416110553.89216507426.3504609207-13657027.02116128646.85816385596.4474609807-13777233.04616146894.38616262536.9534610407-13896508.54616165287.01616138256.7624611007-14014845.33316183815.20216012764.8664611607-14132235.35116202469.31815886070.3514612207-14248670.67316221239.65815758182.3944612807-14364143.50516240116.43715629110.2674613407-14478646.18616259089.79615498863.3364614007-14592171.19116278149.79915367451.0574614607-14704711.12616297286.43915234882.9784615207-14816258.73916316489.63615101168.7384615807-14926806.91316335749.24214966318.0674616407-15036348.66816355055.04214830340.7834617007-15144877.16616374396.75514693246.7954617607-15252385.71016393764.03614555046.0994618207-15358867.74216413146.47614415748.7814618807-15464316.85116432533.61014275365.0134619407-15568726.76516451914.91114133905.0534620007-15672091.35816471279.79713991379.2457心得体会与建议此次作业是本学期“GPS卫星导航定位技术与方法”课程的第一次编程实践作业。在学习了教材第三章“卫星轨道运动及卫星坐标计算”之后，初步对卫星坐标计算有一定认识，担不是很深刻。开始此次作业时，又温习了一遍这个部分的内容，进一步理解勒表示卫星位置的一些参数和量，还有计算的过程和步骤，然后根据已知条件从GPS导航文件中提取出计算所需的广播星历计算参数，按照程序导入文件要求的相应格式输入到txt文本中，作为导入已知数据准备。在编程的过程中，由于去年摄影测量课程的编程基础，感觉编写的代码还算得心应手，没有遇到什么特别大的困难。不过有些小地方出错，例如：1.忘记调用导入已知数-10-据的子程序；2.还有计算过程中，求弧度值的反正切涉及的函数选择问题，一开始使用的是tan（）函数，但是这个函数不能很好解决弧度区间的问题，以至计算错误，后来改用了tan2（