武汉大学卫星导航原理课程作业——广播星历与卫星位置盘算 ...

武汉大学卫星导航原理课程作业
采用C++语言编译，实现rinex3.04格式的广播星历文件读取和卫星的位置解算。

• 作业内容
  

• 根据所给定的广播星历文件，盘算当天多颗GPS卫星（3颗及以上）在24h内每隔30s在瞬时地心地固系下的坐标。
  
• 将上述盘算得到的GPS卫星位置与精密星历对比（15min间隔），评估广播星历盘算出的卫星位置的误差，绘制出卫星的轨道误差时序图。
  

• 给定文件说明
  

• 广播星历文件：BRDC00IGS_R_20213480000_01D_MN.rnx
  
• 精密星历文件：WUM0MGXFIN_20213480000_01D_15M_ORB.SP3
  
• rinex 04格式说明文档。给出的广播星历文件格式为rinex 04，参数含义参照说明文档附录A5、A6。
  

一、广播星历文件的读取

1.首先读取所有内容，存储到lines容器里面，并记录下头文件竣事的位置，便于后续获取数据

1. vector<string> lines;
2. string end_of_header = "                                                            END OF HEADER";
3. int end_of_line = 0;//文件头结束的位置
4. void readfile(string filename)
5. {
6.         int i = 0;
7.         ifstream infile(filename);
8.         if (!infile.is_open())
9.         {
10.                 cerr << "open error!" << endl;
11.         }
12.         string line="0";
13.         while (getline(infile, line))
14.         {
15.                 if (line == end_of_header)
16.                 {
17.                         end_of_line = i;
18.                 }
19.                 lines.push_back(line);
20.                 i = i + 1;
21.         }
22.         infile.close();
23. }//读取广播星历文件

复制代码

2.界说结构体并进行数据分配

1. class Data
2. {
3. private:
4.         struct data_body
5.         {
6.                 string sPRN;//卫星号
7.                 int TOC_Y;
8.                 int TOC_M;
9.                 int TOC_D;
10.                 int TOC_H;
11.                 int TOC_Min;
12.                 int TOC_Sec;//卫星钟的参考时刻，年月日...
13.                 double sa0;//钟差
14.                 double sa1;//钟速
15.                 double sa2;//钟漂
17.                 //广播轨道-1
18.                 double IODE;//星历发布时间
19.                 double Crs;
20.                 double deltan;
21.                 double M0;
23.                 //广播轨道-2
24.                 double Cuc;
25.                 double e;//轨道偏心角
26.                 double Cus;
27.                 double sqrtA;//长半轴平方根
29.                 //广播轨道-3
30.                 double TOE;//星历的参考时刻（GPS周内秒）
31.                 double Cic;
32.                 double OMEGA;//升交点赤经
33.                 double Cis;
35.                 //广播轨道-4
36.                 double i0;//参考时刻的轨道倾角
37.                 double Crc;
38.                 double omega;//近地点角距
39.                 double deltaomega;//
41.                 //广播轨道-5
42.                 double IDOT;
43.                 double L2code;
44.                 double GPSweek;//GPS周数
45.                 double L2flag;
47.                 //广播轨道-6
48.                 double aACC;//卫星精度
49.                 double sHEA;//卫星健康状态
50.                 double TGD;//TGD
51.                 double IODC;//钟的数据龄期
53.                 //广播轨道-7
54.                 double TTN;//电文发送时刻
55.                 double fit;//拟合区间
56.                 double none1;
57.                 double none2;//备用
58.         };
59.         struct position
60.         {
61.                 int GPSweek;
62.                 double GPSsec;
63.                 double x;
64.                 double y;
65.                 double z;//地心地固坐标系下的坐标
66.         };
67. public:
68.         vector <data_body> data;
69.         vector <position> data_P;
70.         void getdata()
71.         {
72.                 for (int j = end_of_line + 1; j <= size(lines) - 1; j = j + 8)
73.                 {
74.                         data_body data1;
75.                         istringstream iss(lines[j]);
76.                         iss >> data1.sPRN;
77.                         iss >> data1.TOC_Y >> data1.TOC_M >> data1.TOC_D >> data1.TOC_H >> data1.TOC_Min >> data1.TOC_Sec;
78.                         iss>> data1.sa0 >> data1.sa1 >> data1.sa2;
79.                         istringstream iss1(lines[j + 1]);
80.                         iss1 >> data1.IODE >> data1.Crs >> data1.deltan >> data1.M0;
81.                         istringstream iss2(lines[j + 2]);
82.                         iss2 >> data1.Cuc >> data1.e >> data1.Cus >> data1.sqrtA;
83.                         istringstream iss3(lines[j + 3]);
84.                         iss3 >> data1.TOE >> data1.Cic >> data1.OMEGA >> data1.Cis;
85.                         istringstream iss4(lines[j + 4]);
86.                         iss4 >> data1.i0 >> data1.Crc >> data1.omega >> data1.deltaomega;
87.                         istringstream iss5(lines[j + 5]);
88.                         iss5 >> data1.IDOT >> data1.L2code >> data1.GPSweek >> data1.L2flag;
89.                         istringstream iss6(lines[j + 6]);
90.                         iss6 >> data1.aACC >> data1.sHEA >> data1.TGD >> data1.IODC;
91.                         istringstream iss7(lines[j + 7]);//逐行解析数据
92.                         vector<double> row;//存储数据的临时容器
93.                         double value=0;
94.                         while (iss7 >> value)
95.                         {
96.                                 row.push_back(value);
97.                         }
98.                         while (row.size() < 4)
99.                         {
100.                                 row.push_back(0);
101.                         }//数据不足4个时，补齐为0
102.                         data1.TTN = row[0];
103.                         data1.fit = row[1];
104.                         data1.none1 = row[2];
105.                         data1.none2 = row[3];
106.                         data.push_back(data1);//将解析完的数据添加到容器里面
107.                 }
108.         }

复制代码

二、卫星位置解算

1.盘算步骤：

  

2.时间转换和最佳历元选取
需要将世界时转换成GPS时，同时需要查找时间差最近的历元，以此历元的数据进行卫星位置解算，代码如下：

1. double get_GPST(int Y, int M, int D, int H, int min, double S)
2. {
3.         double y = 0;
4.         double m = 0;
5.         if (M <= 2)
6.         {
7.                 y = Y - 1;
8.                 m = M + 12;
9.         }
10.         else
11.         {
12.                 y = Y;
13.                 m = M;
14.         }
15.         double UT = 0;//世界时
16.         double MJD = 0;//简化儒略日
17.         UT = H + min / 60.0 + S / 3600.0;
18.         MJD = int(365.25*y) + int(30.6001*(m + 1)) + D + UT / 24 - 679019;
19.         int GPSweek = 0;
20.         double secofweek = 0;
21.         GPSweek = int((MJD - 44244) / 7);//GPS周
22.         secofweek = (MJD - 44244 - GPSweek * 7) * 86400;//GPS周秒
23.         return secofweek;
24. }//计算GPS周内秒的函数
25. int select_epoch(double secofweek,string PRN)
26.         {
27.                 int n = 0;
28.                 double min_s = 10000;
29.                 double Min;
30.                 for (int i = 0; i <= size(data) - 1; i++)
31.                 {
32.                         if (data[i].sPRN == PRN)
33.                         {
34.                                 Min = fabs(secofweek - data[i].TOE);
35.                                 if (Min <= min_s)
36.                                 {
37.                                         n = i;
38.                                         min_s = Min;
39.                                 }
40.                         }
41.                 }
42.                 return n;
43.         }//选择时间差最小的历元计算
44.         void caculate(int k,double secofweek)
45.         {
46.                 double A = pow(data[k].sqrtA,2);
47.                 double n0 = sqrt(u / pow(A, 3));
48.                 double n = n0 + data[k].deltan;//计算卫星的平均角速度n
49.                 double tk = secofweek - data[k].TOE;//计算规划时间
50.                 if (tk > 302400)
51.                 {
52.                         tk = tk - 604800;
53.                 }
54.                 if (tk < -302400)
55.                 {
56.                         tk = tk + 604800;
57.                 }//考虑跨周的影响
58.                 double Mk = data[k].M0 + n * tk;//计算平近点角
59.                 double E0 = 0;
60.                 double E1 = Mk;
61.                 while (abs(E1-E0)>1e-6)
62.                 {
63.                         E0 = E1;
64.                         E1 = Mk + data[k].e*sin(E0);
65.                 }//计算偏近点角
66.                 double cosv = (cos(E1) - data[k].e) / (1 - data[k].e * cos(E1));
67.                 double sinv = sqrt(1 - pow(data[k].e, 2))*sin(E1) / (1 - data[k].e * cos(E1));
68.                 double Vk = atan2(sinv, cosv);
69.                 double w = Vk + data[k].omega;
70.                 double Ou = data[k].Cuc*cos(2 * w) + data[k].Cus*sin(2 * w);
71.                 double Or = data[k].Crc*cos(2 * w) + data[k].Crs*sin(2 * w);
72.                 double Oi = data[k].Cic*cos(2 * w) + data[k].Cis*sin(2 * w);
73.                 double uk = w + Ou;
74.                 double rk = A * (1 - data[k].e*cos(E1)) + Or;
75.                 double ik = data[k].i0 + Oi + data[k].IDOT*tk;
76.                 double xk1 = rk * cos(uk);
77.                 double yk1 = rk * sin(uk);
78.                 double Ok = data[k].OMEGA + (data[k].deltaomega - Oe)*tk - Oe * data[k].TOE;
79.                 double xk = xk1 * cos(Ok) - yk1 * cos(ik)*sin(Ok);
80.                 double yk = xk1 * sin(Ok) + yk1 * cos(ik)*cos(Ok);
81.                 double zk = yk1 * sin(ik);
82.                 position P;
83.                 P.GPSweek = data[k].GPSweek;
84.                 P.GPSsec = secofweek;
85.                 P.x = xk;
86.                 P.y = yk;
87.                 P.z = zk;
88.                 data_P.push_back(P);//把坐标数据添加到容器里面
89.         }

复制代码

三、结果文件保存

将盘算所得结果按照指定要求保存在文本文件中：

1. void savefile()
2.         {
3.                 ofstream outfile("广播星历.txt");
4.                 outfile <<fixed<< setprecision(3);
5.                 for (int i = 0; i <= size(data_P) - 1; i++)
6.                 {
7.                         outfile << data_P[i].GPSweek << " " << data_P[i].GPSsec << " " << data_P[i].x << " " << data_P[i].y << " " << data_P[i].z << endl;
8.                 }
9.         }//将广播星历的数据写入文件

复制代码

  武汉大学测绘学院在读本科生一枚，本次的内容就分享到这里了，渴望各人多多关注，后面会连续给各人分享更多的优质内容。本次作业的所有测试用例和源代码都可以免费下载，渴望各人点点关注支持一下。
【免费】武汉大学卫星导航原理课程作业-广播星历与卫星位置盘算.zip资源-CSDN文库



免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。