该文章提供MATLAB函数，用于处理GEO和MEO/IGSO卫星数据，包括计算卫星运行参数、轨道坐标等。用户可配置时间参数，函数从TXT文件中提取北斗广播星历数据进行计算。
MATLAB文件主函数
GitHub源程序地址：
https://github.com/brucema49/-.git
GEO卫星数据处理函数
自行配置是否需要year，month等时间参数
function [Xk,Yk,Zk,xk,yk,tk,sr] = GEO(year,month,day,hour,minute,second,a0,a1,a2,IodeIssueofData,Crs,deta_n,M0,Cuc,...
e,Cus,a_sqrt,t0e,Cic,W0,Cis,i0,Crc,w,W_,I_)
%1.计算卫星运行的平均角速度n
u=3.986004418*10^14;
n0=sqrt(u)/(a_sqrt)^3;
n=n0+deta_n;
%2.计算归化时间tk
s1=floor((1461*(year+4800+ floor((month-14)/12 )))/4 );
s2= floor(367*(month-2-floor((month-14)/12 )*12)/12 );
s3=floor(3*floor((year+4900+floor((month-14)/12) )/100)/4 );
jd=day-32075+s1+s2-s3-0.5+hour/24+minute/1440+second/86400;
z=floor((jd-2453736.5)/7);
t_=((jd-2453736.5)/7-z)*604800;
t0c=t_;
deta_t=a0+a1*(t_-t0c)+a2*(t_-t0c);
t=t_-deta_t;
tk=t-t0e;
%3.观测时刻卫星平近点角Mk的计算
Mk=M0+n*tk;
%4.计算偏近点角Ek
E=1;
Ek=Mk;
while E>10^(-6)
E=Ek;
Ek=Mk+e*sin(Ek);
E=abs(Ek-E);
end;
%5.真近点角Vk的计算
Vk=atan((sqrt(1-e*e)*sin(Ek))/(cos(Ek)-e));
%6.升交距角Fk的计算
Fk=Vk+w;
%7.摄动改正项su,sr,si的计算
su=Cuc*cos(2*Fk)+Cus*sin(2*Fk);
sr=Crc*cos(2*Fk)+Crs*sin(2*Fk);
si=Cic*cos(2*Fk)+Cis*sin(2*Fk);
%8.计算经过改正的升交距角uk、卫星矢径rk和轨道倾角ik
uk=Fk+su;
rk=a_sqrt*a_sqrt*(1-e*cos(Ek))+sr;
ik=i0+si+I_*tk;
%9.计算卫星在轨道平面直角坐标系的坐标
xk=rk*cos(uk);
yk=rk*sin(uk);
%10计算观测时刻升交点经度Wk
we=7.2921150*10^(-5);
Wk=W0+W_*tk-we*t0e;
%11计算GEO卫星在自定义坐标系中的空间直角坐标
XGk=xk*cos(Wk)-yk*cos(ik)*sin(Wk);
YGk=xk*sin(Wk)+yk*cos(ik)*cos(Wk);
ZGk=yk*sin(ik);
%12计算GEO卫星在CGCS2000地固坐标系中的空间直角坐标系
Rx=[1 0 0;0 cosd(-5) sind(-5);0 -sind(-5) cosd(-5)];
Rz=[cos(we*tk) sin(we*tk) 0;-sin(we*tk) cos(we*tk) 0;0 0 1 ];
ju=Rz*Rx*[XGk;YGk;ZGk];
Xk=ju(1,1);
Yk=ju(2,1);
Zk=ju(3,1);
%禁止一切商业转载，网络转载请注明出处
%https://blog.csdn.net/yitian_00/category_12244341.html?spm=1001.2014.3MEO/IGSO卫星数据处理函数
自行配置是否需要year，month等时间参数
function [Xk,Yk,Zk,xk,yk,tk] = MEO_IGSO(year,month,day,hour,minute,second,a0,a1,a2,IodeIssueofData,Crs,deta_n,M0,Cuc,...
e,Cus,a_sqrt,t0e,Cic,W0,Cis,i0,Crc,w,W_,I_)
%1.计算卫星运行的平均角速度n
u=3.986004418*10^14;
n0=sqrt(u)/(a_sqrt)^3;
n=n0+deta_n;
%2.计算归化时间tk
s1=floor((1461*(year+4800+ floor((month-14)/12 )))/4 );
s2= floor(367*(month-2-floor((month-14)/12 )*12)/12 );
s3=floor(3*floor((year+4900+floor((month-14)/12) )/100)/4 );
jd=day-32075+s1+s2-s3-0.5+hour/24+minute/1440+second/86400;
z=floor((jd-2453736.5)/7);
t_=((jd-2453736.5)/7-z)*604800;
t0c=t_;
deta_t=a0+a1*(t_-t0c)+a2*(t_-t0c);
t=t_-deta_t;
tk=t-t0e;
%3.观测时刻卫星平近点角Mk的计算
Mk=M0+n*tk;
%4.计算偏近点角Ek
E=1;
Ek=Mk;
while E>10^(-6)
E=Ek;
Ek=Mk+e*sin(Ek);
E=abs(Ek-E);
end;
%5.真近点角Vk的计算
Vk=atan((sqrt(1-e*e)*sin(Ek))/(cos(Ek)-e));
%6.升交距角Fk的计算
Fk=Vk+w;
%7.摄动改正项su,sr,si的计算
su=Cuc*cos(2*Fk)+Cus*sin(2*Fk);
sr=Crc*cos(2*Fk)+Crs*sin(2*Fk);
si=Cic*cos(2*Fk)+Cis*sin(2*Fk);
%8.计算经过改正的升交距角uk、卫星矢径rk和轨道倾角ik
uk=Fk+su;
rk=a_sqrt*a_sqrt*(1-e*cos(Ek))+sr;
ik=i0+si+I_*tk;
%9.计算卫星在轨道平面直角坐标系的坐标
xk=rk*cos(uk);
yk=rk*sin(uk);
%10计算观测时刻升交点经度Wk
we=7.292115*10^(-5);
Wk=W0+(W_-we)*tk-we*t0e;
%11
Xk=xk*cos(Wk)-yk*cos(ik)*sin(Wk);
Yk=xk*sin(Wk)+yk*cos(ik)*cos(Wk);
Zk=yk*sin(ik);
%禁止一切商业转载，网络转载请注明出处
%https://blog.csdn.net/yitian_00/category_12244341.html?spm=1001.2014.3数据文件准备工作
请将你的北斗广播星历提取到TXT文件，并确保数据块从第1行开始（如若不为第1行，自行变更主文件line参数）
主函数数据提取文件
文件提取部分参考：MATLAB--读取广播星历的导航文件_导航电文 matlab_CUG扛把子的博客-CSDN博客
%禁止一切商业转载，网络转载请注明出处
%https://blog.csdn.net/yitian_00/category_12244341.html?spm=1001.2014.3
%打开文件
BDS=fopen('你的文件路径','r');
%判断文件是否读取成功
if BDS<0
error('文件读取失败。');
end
%文件行数
line=0;
%参数过多，在行末以...表示连接下一行
satellite=struct('PRN','Year','Month','Day','Hour','Minute',...
'Second','ClockBias','ClockDirft','ClockDirftRate','IodeIssueofData',...
'Crs','Delta_n','M0','Cuc','e','Cus' ,'sqrt_A',...
'toe','Cic','OMEGA','Cis','i0','Crc','omega', 'OMEGADOT','IDOT',...
'CodesOnL2Channel','GPSWeek','L2_p_data_flag','SV_accuracy','SV_health',...
'TGD','IODC_Issue_of_Data','Transmission_time_of_message','Fit_interval');
satellite=repmat(satellite,[1000 40]);
%读取数据
i=1;%卫星数量
g=1;m=1;tk=zeros(200,1);
Mx=zeros(200,1);My=zeros(200,1);Mz=zeros(200,1);
Gx=zeros(50,1);Gy=zeros(50,1);Gz=zeros(50,1);
while ~feof(BDS)
str=fgetl(BDS);
str=strrep(str,'D','e');
line=line+1;
%直接使用动态结构体satellite
if mod(line,8)==1
satellite(i).PRN=str2double(str(2:3));
satellite(i).Year=str2double(str(5:8));
satellite(i).Month=str2double(str(10:11));
satellite(i).Day=str2double(str(13:14));
satellite(i).Hour=str2double(str(16:17));
satellite(i).Minute=str2double(str(19:20));
satellite(i).Second=str2double(strrep(str(22:23),' ',''));
satellite(i).ClockBias=str2double(strrep(str(24:42),' ',''));
satellite(i).ClockDirft=str2double(strrep(str(43:61),' ',''));
satellite(i).ClockDirftRate=str2double(strrep(str(62:80),' ',''));
end
if mod(line,8)==2
satellite(i).IodeIssueofData=str2double(strrep(str(5:23),' ',''));
satellite(i).Crs=str2double(strrep(str(24:42),' ',''));
satellite(i).Delta_n=str2double(strrep(str(43:61),' ',''));
satellite(i).M0=str2double(strrep(str(62:80),' ',''));
end
if mod(line,8)==3
satellite(i).Cuc=str2double(strrep(str(5:23),' ',''));
satellite(i).e=str2double(strrep(str(24:42),' ',''));
satellite(i).Cus=str2double(strrep(str(43:61),' ',''));
satellite(i).sqrt_A=str2double(strrep(str(62:80),' ',''));
end
if mod(line,8)==4
satellite(i).toe=str2double(strrep(str(5:23),' ',''));
satellite(i).Cic=str2double(strrep(str(24:42),' ',''));
satellite(i).OMEGA=str2double(strrep(str(43:61),' ',''));
satellite(i).Cis=str2double(strrep(str(62:80),' ',''));
end
if mod(line,8)==5
satellite(i).i0=str2double(strrep(str(5:23),' ',''));
satellite(i).Crc=str2double(strrep(str(24:42),' ',''));
satellite(i).omega=str2double(strrep(str(43:61),' ',''));
satellite(i).OMEGADOT=str2double(strrep(str(62:80),' ',''));
end
if mod(line,8)==6
satellite(i).IDOT=str2double(strrep(str(5:23),' ',''));
satellite(i).CodesOnL2Channel=str2double(strrep(str(24:42),' ',''));
satellite(i).GPSWeek=str2double(strrep(str(43:61),' ',''));
satellite(i).L2_p_data_flag=str2double(strrep(str(62:80),' ',''));
end
if mod(line,8)==7
satellite(i).SV_accuracy=str2double(strrep(str(5:23),' ',''));
satellite(i).SV_health=str2double(strrep(str(24:42),' ',''));
satellite(i).TGD=str2double(strrep(str(43:61),' ',''));
satellite(i).IODC_Issue_of_Data=str2double(strrep(str(62:80),' ',''));
end
if mod(line,8)==0
satellite(i).Transmission_time_of_message=str2double(strrep(str(5:23),' ',''));
satellite(i).Fit_interval=str2double(strrep(str(24:42),' ',''));
%一个数据块已读完,依据PRN对卫星轨道算法分类
if satellite(i).PRN==1||satellite(i).PRN==2||satellite(i).PRN==3||satellite(i).PRN==4||satellite(i).PRN==5||satellite(i).PRN==59||...
satellite(i).PRN==60||satellite(i).PRN==61
[Xk(i),Yk(i),Zk(i),xk(i),yk(i),tk(i),sr(i)]=GEO( satellite(i).Year,satellite(i).Month,satellite(i).Day,satellite(i).Hour,satellite(i).Minute ,...
satellite(i).Second,satellite(i).ClockBias,satellite(i).ClockDirft,satellite(i).ClockDirftRate, ...
satellite(i).IodeIssueofData,satellite(i).Crs , satellite(i).Delta_n, satellite(i).M0, ...
satellite(i).Cuc, satellite(i).e, satellite(i).Cus, satellite(i).sqrt_A, ...
satellite(i).toe, satellite(i).Cic, satellite(i).OMEGA, satellite(i).Cis, ...
satellite(i).i0, satellite(i).Crc, satellite(i).omega,satellite(i).OMEGADOT,satellite(i).IDOT);
%卫星定轨只需要传到IDOT
Gx(g,1)=Xk(i);Gy(g,1)=Yk(i);Gz(g,1)=Zk(i);g=g+1;
else
[Xk(i),Yk(i),Zk(i),xk(i),yk(i),tk(i)]=MEO_IGSO(satellite(i).Year,satellite(i).Month,satellite(i).Day,satellite(i).Hour,satellite(i).Minute ,...
satellite(i).Second,satellite(i).ClockBias,satellite(i).ClockDirft,satellite(i).ClockDirftRate, ...
satellite(i).IodeIssueofData,satellite(i).Crs , satellite(i).Delta_n, satellite(i).M0, ...
satellite(i).Cuc, satellite(i).e, satellite(i).Cus, satellite(i).sqrt_A, ...
satellite(i).toe, satellite(i).Cic, satellite(i).OMEGA, satellite(i).Cis, ...
satellite(i).i0, satellite(i).Crc, satellite(i).omega,satellite(i).OMEGADOT,satellite(i).IDOT);
Mx(m,1)=Xk(i);My(m,1)=Yk(i);Mz(m,1)=Zk(i);m=m+1;
end
position(i,1)=Xk(i);position(i,2)=Yk(i); position(i,3)=Zk(i);
i=i+1;
end
end
disp('位置矩阵为');
disp(position);
%禁止一切商业转载，网络转载请注明出处
%https://blog.csdn.net/yitian_00/category_12244341.html?spm=1001.2014.3禁止一切商业转载，网络转载请注明出处
https://blog.csdn.net/yitian_00/category_12244341.html?spm=1001.2014.3
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
