利用星历数据解算北斗卫星位置

网上已经有了比较多关于如何利用播发的广播星历来解算卫星位置的blog，此文章的目的是为了记录一下学习成果的同时回馈一下。（毕竟看过很多blog了，但是从来没写过）,我在文中用到的数据是我处理过的，只保留了北斗的数据，如果直接用rinex格式下的星历文件，同样可以按照相同的步骤进行求解，不过在数据处理方面有所不同。

在实现之前首先得下载好广播星历文件并且导入

下载的网址就比较多了，比如说武汉IGS数据中心，广播星历和精密星历都能下载。

下面是具体实现

1.首先对参数初始化

代码目前没有做交互，也没有采用循环（这些都比较容易在后期实现）所以在一开始只能选定一个卫星号prn，并且只能人为修改prn号来更换卫星的选择。在观测时间的选择上，我所读取的广播星历数据的起始时间为BDT 783周259200s，所以我选取的初始观测时间为259200s。计算卫星位置的数量为可以自己选择，由于我下载的精密星历同一卫星仅有95个数据，所以我设定为95个。卫星的时间间隔我选择的是900s，也即15mins。(与精密星历播发数据相匹配方便精度检验)。

clc;
close;
clear;

%1. 读取下载好的广播星历数据
Data_bro=importdata("Nfile-new.csv");

%2.参数初始化
prn=07;             %卫星号 
t=259200;           %观测起始时间
number=95;          %设定计算位置点的数量(当gap=900时，number应小于95，精密星历中只有95条数据)
gap=900;            %时间上的间隔(单位为s)
flag=0;             %选择是否播放动画(flag=0时不播放，flag=1时播放)

2.计算卫星位置

计算卫星位置之前首先要对星历数据进行选择，选择星历要进行两步处理，第一是在星历文件中找到与所设定卫星号prn相同的星历数据，第二步是在第一步的结果中选取当前观测时刻内依旧有效的精度最高的卫星星历数据，由于我们是非实时解算，所以我们可以选取与观测时刻时间间隔最小的星历数据作为我们解算位置所用的数据。也就是说最终选取得到的结果是针对单一卫星在某时刻播发的一条星历数据。

%3.计算卫星位置
for i=1:number      
    t=t+gap;    %每隔gap观测一次
    broadEph_tar = select_eph(t,prn,Data_bro.data);   %选星历
    [Xk(i),Yk(i),Zk(i)] = cal_coordinate(t,prn,broadEph_tar);
end

其中select_eph函数为我设立的选取星历的函数，cal_coordinate为我计算卫星位置的函数。(其实可以把两个函数二合一)其代码如下
select_eph函数:

function broadEph_tar = select_eph(t,prn,Broad_eph)

%1.选择prn相对应的所有星历数据
result_pre=[];
num=1;
for k=1:length(Broad_eph)
    if Broad_eph(k,1)==prn     
        result_pre(num,1)=num;
        result_pre(num,2)=k;
        num=num+1;
    end
end

%2.选择与观测时间相近且已经播发的星历数据  
delta_t=[];
%从对应的广播星历中找到与观测时间最相近的星历数据
for k=1:size(result_pre,1)
    delta_t(k)=t-Broad_eph(result_pre(k,2),3);
    %使用在有效期内的星历数据，将尚未播发的星历数据剔除
    if delta_t(k)<0
        delta_t(k)=1e8;
    end
end
%找到最近播发的星历文件及在数据中所排位置eph_tar,返回所找到的星历文件
[Min,c]=min(delta_t);
eph_tar=result_pre(c,2);
broadEph_tar=Broad_eph(eph_tar,:);
end

cal_coordinate函数：（由于每个人对数据的处理可能不一样，我在下面的代码中省去了我将数据由data_broadEph(也就是所选星历)中读取出来的部分），函数中主要给出的是具体的计算步骤，其实也比较简单。其中有一个对卫星状态的判断，上官网直接看中国卫星导航系统管理办公室测试评估研究中心

function [Xk,Yk,Zk] = cal_coordinate(t,prn,data_broadEph)
%1.解算所用常数
miu=3.986004418e14;
Omega_Dote=7.2921150e-5;

%2.计算轨道半长轴A及卫星平均角速度n0
A=sqrt_A^2;
n0=sqrt(miu/A^3);

%3.计算时间差
% toc=toc_week*3600*24*7+toc_second;
% toe=toe_week*3600*24*7+toe_second;
toe=toe_second;
tk=t-toe;
% if tk>302400
%     tk=tk-604800;
% elseif tk<-302400
%     tk=tk+604800;
% end

%4.改正平均角速度n0
n=n0+delta_n;

%5.计算平近点角
Mk=M0+n*tk;

%6.迭代计算偏近点角
%对开普勒方程进行迭代得到三阶近似式
E0=Mk+e*sin(Mk)+e^2*sin(Mk)*cos(Mk)+0.5*e^3*sin(Mk)*(3*(cos(Mk))^2-1);
% E0=Mk;
E1=E0+e*sin(E0);
while abs(E1-E0)>1e-10
    E0=E1;
    E1=Mk+e*sin(E0);
    delta=abs(E1-E0);
end
Ek=E1;

%7.计算纬度幅角
% vk_sin=sqrt(1-e^2)*sin(Ek)/(1-e*cos(Ek));
% vk_cos=cos(Ek)-e/(1-e*cos(Ek));
vk=2*atan2(sqrt((1+e)/(1-e))*tan(Ek/2),1);
phi_k=vk+omega;

%8.计算纬度幅角、径向、轨道倾角改正项
delta_uk=Cus*sin(2*phi_k)+Cuc*cos(2*phi_k);
delta_rk=Crs*sin(2*phi_k)+Crc*cos(2*phi_k);
delta_ik=Cis*sin(2*phi_k)+Cic*cos(2*phi_k);

%9.改正纬度幅角、径向、轨道倾角
uk=phi_k+delta_uk;
rk=A*(1-e*cos(Ek))+delta_rk;
ik=i0+Idot*tk+delta_ik;

%10.计算卫星再轨道平面内的坐标
xk=rk*cos(uk);
yk=rk*sin(uk);

%% 计算卫星在BDCS坐标系下的坐标
% 1.首先对卫星状态进行判断
%如果卫星为为非正常状态(在轨实验/在轨测试)则不返回卫星位置
prn_abnormal=[31,56,57,18];
prn_Geo=[1,2,3,4,5,59,69];
% if ismember(prn,prn_abnormal)==1
%     return
if ismember(prn,prn_Geo)==1   %GEO
    %(1)计算历元升交点经度(惯性系)
    Omega_k=Omega0+OmegaDot*tk-Omega_Dote*toe;
    %(2)计算GEO卫星在自定义坐标系中的坐标
    XGk=xk*cos(Omega_k)-yk*cos(ik)*sin(Omega_k);
    YGk=xk*sin(Omega_k)+yk*cos(ik)*cos(Omega_k);
    ZGk=yk*sin(ik);
    %(3)计算GEO卫星在BDCS中的坐标
    Xk=cos(Omega_Dote*tk)*XGk+sin(Omega_Dote*tk)*cos(-5/180*pi)*YGk+sin(Omega_Dote*tk)*sin(-5/180*pi)*ZGk;
    Yk=-sin(Omega_Dote*tk)*XGk+cos(Omega_Dote*tk)*cos(-5/180*pi)*YGk+cos(Omega_Dote*tk)*sin(-5/180*pi)*ZGk;
    Zk=-sin(-5/180*pi)*YGk+cos(-5/180*pi)*ZGk;
else
    %计算MEO/IGSO卫星坐标
    Omega_k=(Omega0+(OmegaDot-Omega_Dote)*tk-Omega_Dote*toe);   %计算升交点经度(地固系)
    Xk=xk*cos(Omega_k)-yk*cos(ik)*sin(Omega_k);
    Yk=xk*sin(Omega_k)+yk*cos(ik)*cos(Omega_k);
    Zk=yk*sin(ik);
end

3.绘制卫星星下点轨迹

首先要将卫星的坐标在大地坐标系下表示，在转换之后要将经纬度由弧度制转换为角度制。其中的绘制动图是为了更好的展示星下点轨迹而增添的一个功能，可以一个一个点逐点将星下点轨迹绘制出来，可以让大家知道星下点的起始点。

%4. 绘制卫星星下点轨迹
%4.1 转换坐标 
pos=[Xk;Yk;Zk]';
% plot3(pos(:,1),pos(:,2),pos(:,3));  grid on
for i=1:length(pos)
    [B(i),L(i),H(i)] = xyz2blh(Xk(i),Yk(i),Zk(i));   %将空间直角坐标转换成大地坐标
end
B=B'*180/pi;
L=L'*180/pi;
h = geoshow('landareas.shp', 'FaceColor', [0.5 0.8 0.3]);  %导入世界地图
hold on;    grid on;
plot(L,B);
%4.2 绘制动图
if flag==1
    at = animatedline('Color','b','Marker','*','MarkerSize',3);
%     avi = VideoWriter('orbit.avi');
%     open(avi);
    Mo=moviein(length(pos));
    for k=1:length(pos)
        addpoints(at,L(k),B(k));
        Mo(k)=getframe;
    end
    drawnow;
    movie(Mo);
%     close(avi);
end

4.读入精密星历计算精度

这一步比较简单，首先对时间间隔进行一个判断，如果间隔为15分钟才会读取精密星历，来计算精度。

if gap==900    %如果时间间隔为900s也即15分钟，才导入精密星历进行误差计算
    Data_pre=importdata("SP3.mat");
    Data_pre=Data_pre(:,1:4);        %只保留卫星坐标
    
    %1.根据所定prn寻找符合的精密星历数据
    if prn<10
        PRN="PC"+"0"+num2str(prn);
    else
        PRN="PC"+num2str(prn);
    end
    k=[];
    num=0;
    for i=1:length(Data_pre)        %找到PRN相匹配的卫星
        if Data_pre{i,1}==PRN
            num=num+1;
            k(num)=i;
        end
    end
    
    PreciseEph=Data_pre(k(2:end),2:4);   %由于未观测起始时间的卫星所以剔除第一个数据
    PreciseEph=cell2mat(PreciseEph);
    error=pos-PreciseEph(1:length(pos),:)*1000;     %精密星历中的坐标数据单位为KM
    error_mean=mean(error);
    
    %2. 绘制图像
    figure;hold on;
    plot(1:length(error),error(:,1));
    plot(1:length(error),error(:,2));
    plot(1:length(error),error(:,3));
    grid on;
    legend("X坐标误差","Y坐标误差","Z坐标误差");
    ylabel("误差大小(m)");
end

由误差图可以看到，由本程序实现的卫星位置计算，精度在米级，误差在5m以内，比较正常，说明上诉解算流程基本正确。
到这里基本上就实现了我们的目的，希望能给大家带来一定的帮助！如果想要通过卫星位置的解算进一步实现用户位置的解算，可以阅读我的另一篇文章