基于广播星历的北斗定位解算原理（基于C语言和MATLAB实现）

一、利用广播星历解算卫星位置

网上已经有很多解算北斗卫星位置的资料和算法，此文章的目的只是记录自己的学习过程，并且积累一些写文章博客的经验。本文先用C语言解算卫星位置，再用MATLAB绘出卫星三维坐标图。本篇博客所使用的资料和文件都是网络上公开发表且可以找到的资料文件。16参数所解算的文件是RINEX3.04标准；我做的18参数模型采用的是国内的格式，故解算的文件是RINEX3.03标准。

以下是实现方法

我所参考的是《2017中国卫星导航系统管理办公室》和《2019中国卫星导航系统管理办公室》文件，链接放在下面，大家感兴趣的可自行下载，本文所引用的公式及截图均来自此文件。

---

2017版：http://www.beidou.gov.cn/yw/xwzx/201712/W020171226734109234149.pdf

2019版：http://www.beidou.gov.cn/xt/gfxz/201902/P020190227592987952674.pdf

---

一、16参数算法

一、前言

16参数解算步骤每次只能解算其中一颗卫星的位置；18参数结算步骤可以直接导入完整的文件，并对其中的北斗卫星进行解算并存入文件中，方便后续导入其他软件做数据研究。

无论是16参数或者是18参数版本，在main函数中参数的顺序是按照（1+6+9）即（时间参数+6个开普勒轨道根数+9个摄动参数）声明的，18参数只是加了两个变化率参数。

卫星型号具体是MEO/GEO/IGSO，可以参考北斗卫星导航系统办公室发表的文件。如下图所示

16/18参数在广播星历中各自所对应的位置在网上均可找到，此处不做赘述。

 

二、实现步骤（具体每一步讲解会在代码中体现并带有注释）

三、C语言源码

//**************************************** Message ***********************************//
//技术交流：folkore0118@163.com
//关注CSDN博主：“彼稷”
//Author: 彼稷 
//All rights reserved	                               
//----------------------------------------------------------------------------------------
// Target Devices: 		DELL-G15
// Tool Versions:       DEVC++5.11
// File name:           Calculation_BDS_Position
// Last modified Date:  2023年10月12日20:00:00
// Last Version:        
// Descriptions:        广播星历&北斗卫星位置 
//----------------------------------------------------------------------------------------
//****************************************************************************************//
#include 
#include 
#include 
#define	Mu			3.986004418e14
#define	Omega_e		7.292115e-5
#define	Pi			3.1415926535898
typedef	struct{
	char	PRN[10];
	int		year;
	int		month;
	int		day;
	int		hour;
	int		minute;
	double	second;
}TIME;
	
	/*--------------------------------------------------------------------------------
	get the BDS time
	*/
int GetsBDSTime(char PRN[10],int year,int month,int day,int hour,int minute,double second,int *WeekNum,double *SecondOfWeek){
	int	DayOfYear	=	0;
	int	DayOfMonth	=	0;
//	/*-------------------------------------
//	year
//	
	for(int DayOfBDS = 2006;DayOfBDS ");
	scanf("%lf",&t);
//	delta_t	=	a0	+	a1 * ((t * 60 + SecondOfWeek) - toc)	+	a2 * pow(((t * 60 + SecondOfWeek) - toc),2);
	trans_t	=	t * 60	+	SecondOfWeek;
	tk	=	trans_t	-	toe	-	14;
	if(tk > 302400){
		tk	=	tk	-	604800;
	}
	else if(tk 撇 	coordinate--->坐标 
		*/
		double	ApostOfX,ApostOfY,ApostOfZ; 
		double	ResultX,ResultY,ResultZ,Resulta0;
		ApostOfX	=	Xk * cos(OMEGA_k)	-	Yk * cos(ik) * sin(OMEGA_k);
		ApostOfY	=	Xk * sin(OMEGA_k)	+	Yk * cos(ik) * cos(OMEGA_k);
		ApostOfZ	=	Yk * sin(ik);
		printf("\tCoordinateOfX	=	%.20lf\n",ApostOfX);
		printf("\tCoordinateOfY	=	%.20lf\n",ApostOfY);
		printf("\tCoordinateOfZ	=	%.20lf\n",ApostOfZ);
		ResultX		=	ApostOfX / 1000.0;
		ResultY		=	ApostOfY / 1000.0;
		ResultZ		=	ApostOfZ / 1000.0;
		Resulta0	=	a0 * 1000000.0;
		printf("\tResultX		=	%.6lf\n",ResultX);
		printf("\tResultY		=	%.6lf\n",ResultY);
		printf("\tResultZ		=	%.6lf\n",ResultZ);
		printf("\tResulta0		=	%.6lf\n",Resulta0);
	/*-----------------------------------------------------------------------------------------------
	将数据写入文件并保存 
	*/
	printf("/-----------------------------------------------------------------------------------------/\n");
	fprintf(pfileWrite,"%s %.6lf %.6lf %.6lf %.6lf\n",time.PRN,ResultX,ResultY,ResultZ,Resulta0);
	printf("/-----------------------------------------------------------------------------------------/\n");
	}
	fclose(pfile);
	fclose(pfileWrite);
	return 0;
}

值得注意的是，toc转到BDT是从2006年1月1日（周日）作为起始点；因为BDT和GPST相差14秒，所以计算归化时间tk时，需要人为减去14秒，调整到和GPST同步。计算tk时，t相当于是对卫星位置做预测，且北斗卫星每隔一个小时更新一次星历，故t的取值小于60分钟。

二、18参数算法

一、前言

18参数模型相对比16参数改动不是特别大，一个是公式的变动，另一个是广播星历对应参数的位置部分有变动。

公式的变动主要体现在下图，具体的还是参考中国卫星导航系统管理办公室提供的官方文件。

二、实现步骤（具体每一步讲解会在代码中体现并带有注释）

三、C语言源码

本代码仅适用于计算北斗卫星国内格式的文件，国外格式的文件需要做相对应的修正。

本代码可以直接读取整个大文件，并计算出所有北斗卫星的位置。

//**************************************** Message ***********************************//
//技术交流：folkore0118@163.com
//关注CSDN博主：“彼稷”
//Author: 彼稷 
//All rights reserved	                               
//----------------------------------------------------------------------------------------
// Target Devices: 		DELL-G15
// Tool Versions:       DEVC++5.11
// File name:           Calculation_BDS_Position
// Last modified Date:  2023年10月18日20:00:00
// Last Version:        
// Descriptions:        广播星历&北斗卫星位置 &18参数 
//----------------------------------------------------------------------------------------
//****************************************************************************************//
#include 
#include 
#include 
#define	Mu			3.986004418e14
#define	Omega_e		7.292115e-5
#define	Pi			3.1415926535898
#define	Aref_MEO	27906100
#define	Aref_GEO	42162200
typedef	struct{
	char	PRN[10];
	int		year;
	int		month;
	int		day;
	int		hour;
	int		minute;
	double	second;
}TIME;
	
	/*--------------------------------------------------------------------------------
	get the BDS time
	*/
int GetsBDSTime(char PRN[10],int year,int month,int day,int hour,int minute,double second,int *WeekNum,double *SecondOfWeek){
	int	DayOfYear	=	0;
	int	DayOfMonth	=	0;
//	/*-------------------------------------
//	year
//	
	for(int DayOfBDS = 2006;DayOfBDS 0\n");
				//scanf("%lf",&t);
				//	delta_t	=	a0	+	a1 * ((t * 60 + SecondOfWeek) - toc)	+	a2 * pow(((t * 60 + SecondOfWeek) - toc),2);
				trans_t	=	t * 60	+	SecondOfWeek;
				tk	=	trans_t	-	toe	-	14;
				if(tk > 302400){
					tk	=	tk	-	604800;
				}
				else if(tk 撇 	coordinate--->坐标 
					*/
					double	ApostOfX,ApostOfY,ApostOfZ; 
					double	ResultX,ResultY,ResultZ,Resulta0;
					ApostOfX	=	Xk * cos(OMEGA_k)	-	Yk * cos(ik) * sin(OMEGA_k);
					ApostOfY	=	Xk * sin(OMEGA_k)	+	Yk * cos(ik) * cos(OMEGA_k);
					ApostOfZ	=	Yk * sin(ik);
					printf("\tCoordinateOfX	=	%.20lf\n",ApostOfX);
					printf("\tCoordinateOfY	=	%.20lf\n",ApostOfY);
					printf("\tCoordinateOfZ	=	%.20lf\n",ApostOfZ);
					ResultX		=	ApostOfX / 1000.0;
					ResultY		=	ApostOfY / 1000.0;
					ResultZ		=	ApostOfZ / 1000.0;
					Resulta0	=	a0 * 1000000.0;
					printf("\tResultX		=	%.6lf\n",ResultX);
					printf("\tResultY		=	%.6lf\n",ResultY);
					printf("\tResultZ		=	%.6lf\n",ResultZ);
					printf("\tResulta0	=	%.6lf\n",Resulta0);
					/*-----------------------------------------------------------------------------------------------
					将数据写入文件并保存 
					*/
					printf("/-----------------------------------------------------------------------------------------/\n");
					fprintf(pfileWrite,"%s %.6lf %.6lf %.6lf %.6lf\n",time.PRN,ResultX,ResultY,ResultZ,Resulta0);
					printf("/-----------------------------------------------------------------------------------------/\n");
				//continue; 
			}
		}
		else{
			continue;
		}
	}
	fclose(pfile);
	fclose(pfileWrite);
	return 0;
}

二、利用MATLAB，并根据18参数解算出来的卫星位置文件绘制三维图

一、代码

MATLAB读文件的位置直接修改为18参数写文件的路径，这样在C程序中直接修改卫星型号，MATLAB就能实时绘图。

%**************************************** Message ***********************************//
%//技术交流：folkore0118@163.com
%//关注CSDN博主：“彼稷”
%//Author: 彼稷
%//All rights reserved	                               
%//----------------------------------------------------------------------------------------
%// Target Devices: 		DELL-G15
%// Tool Versions:       MATLAB R2021a
%// File name:           draw
%// Last modified Date:  2023年10月18日14:00:00
%// Last Version:        
%// Descriptions:        卫星位置&三维图
%//----------------------------------------------------------------------------------------
% 文件路径
filename = '18参数写文件的位置';
% 使用 importdata 函数读取文件
data = importdata(filename, ' ');
% 获取名称、X坐标、Y坐标和Z坐标
name = data.textdata;
x = data.data(:, 1);
y = data.data(:, 2);
z = data.data(:, 3);
% 绘制三维图
figure;
scatter3(x, y, z);
hold on;
% 添加球面
center = [0, 0, 0]; % 球心坐标
radius = 6371; % 球半径
[X, Y, Z] = sphere;
X = radius * X;
Y = radius * Y;
Z = radius * Z ;
surf(X, Y, Z, 'FaceColor', 'none', 'EdgeColor', 'red');
axis equal;
grid on;
xlabel('X');
ylabel('Y');
zlabel('Z');
title('SUM坐标图');

二、三维图

外围蓝色空心的是北斗卫星，中间红色的球面是地球。

根据图中卫星位置来看，在南极和北极上空均没有卫星分布；且IGSO卫星呈“8”字绕行，符合客观事实。

MEO卫星分布：

IGSO卫星分布：

GEO卫星分布：

所有北斗卫星分布：

 

三、最后

本文只是记录学习过程中的一些感想并复盘每阶段的任务，代码仍有优化的空间、内容有谬误之处请大家多多指正。