卫星导航原理课程作业——广播星历与卫星位置计算

武汉大学卫星导航原理课程作业

采用C++语言编译，实现rinex3.04格式的广播星历文件读取和卫星的位置解算。

• 作业内容

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

• 给定文件说明

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

一、广播星历文件的读取

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

vector<string> lines;
string end_of_header = "                                                            END OF HEADER";
int end_of_line = 0;//文件头结束的位置
void readfile(string filename)
{
	int i = 0;
	ifstream infile(filename);
	if (!infile.is_open())
	{
		cerr << "open error!" << endl;
	}
	string line="0";
	while (getline(infile, line))
	{
		if (line == end_of_header)
		{
			end_of_line = i;
		}
		lines.push_back(line);
		i = i + 1;
	}
	infile.close();
}//读取广播星历文件

2.定义结构体并进行数据分配

class Data
{
private:
	struct data_body
	{
		string sPRN;//卫星号
		int TOC_Y;
		int TOC_M;
		int TOC_D;
		int TOC_H;
		int TOC_Min;
		int TOC_Sec;//卫星钟的参考时刻，年月日...
		double sa0;//钟差
		double sa1;//钟速
		double sa2;//钟漂

		//广播轨道-1
		double IODE;//星历发布时间
		double Crs;
		double deltan;
		double M0;

		//广播轨道-2
		double Cuc;
		double e;//轨道偏心角
		double Cus;
		double sqrtA;//长半轴平方根

		//广播轨道-3
		double TOE;//星历的参考时刻（GPS周内秒）
		double Cic;
		double OMEGA;//升交点赤经
		double Cis;

		//广播轨道-4
		double i0;//参考时刻的轨道倾角
		double Crc;
		double omega;//近地点角距
		double deltaomega;//

		//广播轨道-5
		double IDOT;
		double L2code;
		double GPSweek;//GPS周数
		double L2flag;

		//广播轨道-6
		double aACC;//卫星精度
		double sHEA;//卫星健康状态
		double TGD;//TGD
		double IODC;//钟的数据龄期

		//广播轨道-7
		double TTN;//电文发送时刻
		double fit;//拟合区间
		double none1;
		double none2;//备用
	};
	struct position
	{
		int GPSweek;
		double GPSsec;
		double x;
		double y;
		double z;//地心地固坐标系下的坐标
	};
public:
	vector <data_body> data;
	vector <position> data_P;
	void getdata()
	{
		for (int j = end_of_line + 1; j <= size(lines) - 1; j = j + 8)
		{
			data_body data1;
			istringstream iss(lines[j]);
			iss >> data1.sPRN;
			iss >> data1.TOC_Y >> data1.TOC_M >> data1.TOC_D >> data1.TOC_H >> data1.TOC_Min >> data1.TOC_Sec;
			iss>> data1.sa0 >> data1.sa1 >> data1.sa2;
			istringstream iss1(lines[j + 1]);
			iss1 >> data1.IODE >> data1.Crs >> data1.deltan >> data1.M0;
			istringstream iss2(lines[j + 2]);
			iss2 >> data1.Cuc >> data1.e >> data1.Cus >> data1.sqrtA;
			istringstream iss3(lines[j + 3]);
			iss3 >> data1.TOE >> data1.Cic >> data1.OMEGA >> data1.Cis;
			istringstream iss4(lines[j + 4]);
			iss4 >> data1.i0 >> data1.Crc >> data1.omega >> data1.deltaomega;
			istringstream iss5(lines[j + 5]);
			iss5 >> data1.IDOT >> data1.L2code >> data1.GPSweek >> data1.L2flag;
			istringstream iss6(lines[j + 6]);
			iss6 >> data1.aACC >> data1.sHEA >> data1.TGD >> data1.IODC;
			istringstream iss7(lines[j + 7]);//逐行解析数据
			vector<double> row;//存储数据的临时容器
			double value=0;
			while (iss7 >> value)
			{
				row.push_back(value);
			}
			while (row.size() < 4)
			{
				row.push_back(0);
			}//数据不足4个时，补齐为0
			data1.TTN = row[0];
			data1.fit = row[1];
			data1.none1 = row[2];
			data1.none2 = row[3];
			data.push_back(data1);//将解析完的数据添加到容器里面
		}
	}

二、卫星位置解算：

1.计算步骤：

  

2.时间转换和最佳历元选取

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

double get_GPST(int Y, int M, int D, int H, int min, double S)
{
	double y = 0;
	double m = 0;
	if (M <= 2)
	{
		y = Y - 1;
		m = M + 12;
	}
	else
	{
		y = Y;
		m = M;
	}
	double UT = 0;//世界时
	double MJD = 0;//简化儒略日
	UT = H + min / 60.0 + S / 3600.0;
	MJD = int(365.25*y) + int(30.6001*(m + 1)) + D + UT / 24 - 679019;
	int GPSweek = 0;
	double secofweek = 0;
	GPSweek = int((MJD - 44244) / 7);//GPS周
	secofweek = (MJD - 44244 - GPSweek * 7) * 86400;//GPS周秒
	return secofweek;
}//计算GPS周内秒的函数
int select_epoch(double secofweek,string PRN)
	{
		int n = 0;
		double min_s = 10000;
		double Min;
		for (int i = 0; i <= size(data) - 1; i++)
		{
			if (data[i].sPRN == PRN)
			{
				Min = fabs(secofweek - data[i].TOE);
				if (Min <= min_s)
				{
					n = i;
					min_s = Min;
				}
			}
		}
		return n;
	}//选择时间差最小的历元计算
	void caculate(int k,double secofweek)
	{
		double A = pow(data[k].sqrtA,2);
		double n0 = sqrt(u / pow(A, 3));
		double n = n0 + data[k].deltan;//计算卫星的平均角速度n
		double tk = secofweek - data[k].TOE;//计算规划时间
		if (tk > 302400)
		{
			tk = tk - 604800;
		}
		if (tk < -302400)
		{
			tk = tk + 604800;
		}//考虑跨周的影响
		double Mk = data[k].M0 + n * tk;//计算平近点角
		double E0 = 0;
		double E1 = Mk;
		while (abs(E1-E0)>1e-6)
		{
			E0 = E1;
			E1 = Mk + data[k].e*sin(E0);
		}//计算偏近点角
		double cosv = (cos(E1) - data[k].e) / (1 - data[k].e * cos(E1));
		double sinv = sqrt(1 - pow(data[k].e, 2))*sin(E1) / (1 - data[k].e * cos(E1));
		double Vk = atan2(sinv, cosv);
		double w = Vk + data[k].omega;
		double Ou = data[k].Cuc*cos(2 * w) + data[k].Cus*sin(2 * w);
		double Or = data[k].Crc*cos(2 * w) + data[k].Crs*sin(2 * w);
		double Oi = data[k].Cic*cos(2 * w) + data[k].Cis*sin(2 * w);
		double uk = w + Ou;
		double rk = A * (1 - data[k].e*cos(E1)) + Or;
		double ik = data[k].i0 + Oi + data[k].IDOT*tk;
		double xk1 = rk * cos(uk);
		double yk1 = rk * sin(uk);
		double Ok = data[k].OMEGA + (data[k].deltaomega - Oe)*tk - Oe * data[k].TOE;
		double xk = xk1 * cos(Ok) - yk1 * cos(ik)*sin(Ok);
		double yk = xk1 * sin(Ok) + yk1 * cos(ik)*cos(Ok);
		double zk = yk1 * sin(ik);
		position P;
		P.GPSweek = data[k].GPSweek;
		P.GPSsec = secofweek;
		P.x = xk;
		P.y = yk;
		P.z = zk;
		data_P.push_back(P);//把坐标数据添加到容器里面
	}

三、结果文件保存

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

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

武汉大学测绘学院在读本科生一枚，本次的内容就分享到这里了，希望大家多多关注，后面会持续给大家分享更多的优质内容