基于MATLAB的BDS卫星轨道仿真

 一、实现功能

本代码实现了以下几个功能

• 根据广播星历信息提取轨道六根数和摄动参数。

• 计算并绘制BDS GEO卫星、IGSO卫星和MEO卫星的无摄轨道平面。

• 计算并绘制MEO卫星在无摄全轨道平面的位置分布。

• 添加摄动参数，推算两小时之内IGSO卫星位置和速度变化，并绘制位置和速度随时间变化曲线。

 二、了解RINEX格式的BDS广播星历格式和内容

以下是一段北斗星历文件

        注意图中END OF HEADER，此字符串以上部分我们不做赘述，可自行查找，我们只需要考虑此字符串以下的部分。

        要想通过卫星星历文件提取数剧，我们就要了解星历格式每个部分的含义，观察可以发现，文档中的数据是一块一块的，观察每一块数据，可得，它是一个8×5的表格，具体每一个数据对应的意义可以看下表：

卫星编号	时间
		轨道半径改正项Crs	平均角速度改正项delta_n	平近点角M0
	升交点角距改正项Cuc	轨道偏心率e	升交点角距改正项Cus	轨道长半轴平方根sqrta
	星历的参考时刻t_oe	轨道倾角改正项Cic	升交点经度omega	轨道倾角改正项Cis
	轨道倾角i0	轨道半径改正项Crc	近地点角距w	升交点赤经变化omega_rate
	轨道倾角的变率i_rate

 三、代码实现

提取星历文件

fileid = fopen("tarc2440.23b");
global sate;

%读取head部分
while 1
    l1ne=fgetl(fileid);
    answer = strfind(l1ne,"END OF HEADER");
    if ~isempty(answer)
        break
    end
end

%保存星历参数数据于sate元胞数组
while 1
    Data_1=textscan(fileid,'C%d 2023 09 01 %d 00 00 %f32 %f32 %f32');
    Data_2=textscan(fileid,'%f32 %f32 %f32 %f32');    
    flag = Data_1{1,1};
    t = Data_1{1,2};
    %1   2      3   4 5  6   7     8         9        10
    %e,sqrta,omega,i0,w,M0,t_oe,delat_n,omega_rate,i_rate,
    % 11  12  13  14  15  16  17
    %crs,cuc,cus,cic,cis,crc,niu
    sate{flag,t+1} = [Data_2{1,2}(2);Data_2{1,4}(2);Data_2{1,3}(3);Data_2{1,1}(4);Data_2{1,3}(4);...
        Data_2{1,4}(1);Data_2{1,1}(3);Data_2{1,3}(1);Data_2{1,4}(4);Data_2{1,1}(5);Data_2{1,2}(1);...
        Data_2{1,1}(2);Data_2{1,3}(2);Data_2{1,2}(3);Data_2{1,4}(3);Data_2{1,2}(4)];
    
    if flag == 62 && t == 23
        break
    end
end
fclose(fileid);

        在此代码中，我们首先将读取文件的指针指到“END OF HEADER”字符串下面的一行，也就是对应的第一块卫星数据，如图所示：

        之后，通过textscan函数匹配一行数据，进而匹配一段数据，提取星历文件的所有参数，之后，将我们需要的轨道六要素和摄动参数等放到sate元胞数组中，并且注释备注sate里对应的数值的物理意义：

         上图为sate变量使用的含义，如：sate{1，3}（5）代表的是1号卫星在第三时刻（2小时）的近地点角距w（欧米噶）。

计算真近点角

        提取星历数据后，我们来计算真近点角，并且将此变量保存到sate元胞数组的第17个数据位。

具体实现代码如下：

%计算真近点角
for j = 1:1:62
    for i = 0:1:23
        if isempty(sate{j,i+1})
            continue;
        end
        M0 = sate{j,i+1}(6);
        e = sate{j,i+1}(1);
        f = @(E) M0 + e*sin(E);
        [E, n] = iteration_method(f, M0, 1e-90, 500);
        niu(i+1) = 2*atan(sqrt((1+e)/(1-e))*tan(E/2));
        sate{j,i+1}(17) = niu(i+1);
     end
end

        此部分使用到了迭代函数iteration_method,此函数代码如下：

iteration_method.m

% 定义迭代函数
%输入参数：f:方程的函数句柄
%          E0: 初始近似解
%          epsilon: 迭代精度要求
%          max_iter: 最大迭代次数
%输出参数：E：迭代值
%          n：迭代次数

function [E, n] = iteration_method(f, E0, epsilon, max_iter)
     E = E0;  % 初始值
     n = 0;   % 迭代次数
     while (n < max_iter)
         E_next = f(E);  % 迭代公式
         if abs(E_next - E) < epsilon % 判断是否达到指定精度
             break;
         end
         E = E_next;
         n = n + 1;
     end
end

加载地球

        此部分代码可直接使用，目的是为了实现美观，此处不再解释。

load_earth.m

function load_earth()
gridNum = 10;
[x, y, z] = sphere(gridNum);
r = 6400000;
x = x*r;y = y*r;z = z*r;
surf(x,y,z);
earth = imread('landOcean.jpg');
% 高配
clouds = imread('cloudCombined.jpg');

hEarth = surface(x ,y, z);
hEarth.FaceColor = 'texturemap';       
hEarth.EdgeColor = 'none';              
hEarth.EdgeAlpha = 0.3;
hEarth.CData = earth;                

view([80 10]);          
axis equal

绘制三种轨道平面

%绘制三种轨道平面
draw_plane(1);

        此部分用到了draw_plane函数，实现代码如下：

draw_plane.m

%输入参数：flag == 0 绘制所有轨道平面
%          flag == 1 绘制三个不同类型轨道平面

function draw_plane(flag)
global sate;
hold on;
%绘制所有轨道平面
if flag == 0
    for i=1:1:62
        if isempty(sate{i,1})
            continue;
        end
        covert_plane(sate{i,1}(2)^2,sate{i,1}(1),sate{i,1}(4),sate{i,1}(3),sate{i,1}(5));%a,e,i0,omega,w
    end
end
    
%绘制三个不同类型轨道平面
if flag == 1
    for i=[1 7 11]
        covert_plane(sate{i,1}(2)^2,sate{i,1}(1),sate{i,1}(4),sate{i,1}(3),sate{i,1}(5));%a,e,i0,omega,w
    end
end

        此函数输入参数flag目的是确定输出类型，flag == 0，绘制所有轨道平面；flag == 1，绘制卫星号为1、7、11的三个轨道平面，三个卫星分别对应三种不同轨道平面。此部分用到了covert_plane函数，代码部分如下：

 covert_plane.m

%绘制卫星运行轨道椭圆
%输入参数：除真近点角外轨道五要素

function covert_plane(a,e,i0,omega,w)
theta = 0:0.001:2*pi;
rho = a*(1-e^2)./(1+e*cos(theta));
[x,y] = pol2cart(theta,rho); %极坐标转为直角坐标
z = zeros(1,length(x));
R3_omega = [cos(omega),sin(omega),0;-sin(omega),cos(omega),0;0,0,1];
R1_i0 = [1,0,0;0,cos(i0),sin(i0);0,-sin(i0),cos(i0)];
R3_w = [cos(w),sin(w),0;-sin(w),cos(w),0;0,0,1];
R = R3_w*R1_i0*R3_omega;
for i = 1:1:6284
    A = R*[x(i);y(i);z(i)];
    x1(i) = A(1);y1(i) = A(2);z1(i) = A(3);
end
plot3(x1,y1,z1,'Color','green');

        此函数输入卫星除真近点角以外的五个轨道要素，从而可以确定一个椭圆，绘制出的此椭圆，即为对应的卫星运行的轨道平面。

        其中，R为转换矩阵，将赤道平面上的椭圆转换到轨道平面。

绘制MEO位置分布

%绘制MEO位置分布
draw_MEO_local(1,1);

        draw_MEO_local代码如下： 

draw_MEO_local.m

%输入参数：flag == 0 绘制所有MEO位置分布
%          flag == 1 绘制某时刻MEO位置分布
%输入参数：t  在flag == 1时，t为某一时刻

function draw_MEO_local(flag,t)
global sate;
%绘制所有MEO位置分布
if flag == 0
    for j = 1:1:62 
        for i = 1:1:24
            if isempty(sate{j,i})
                continue;
            end
            if(sate_class(sate{j,i}(2),sate{j,i}(4)) == 2)
%                 if j == 11
%                     covert_point(sate{j,i}(2)^2,sate{j,i}(1),sate{j,i}(4),sate{j,i}(3),sate{j,i}(5),sate{j,i}(17),'blue');
%                 else
                        covert_point(sate{j,i}(2)^2,sate{j,i}(1),sate{j,i}(4),sate{j,i}(3),sate{j,i}(5),sate{j,i}(17),'red');
                        hold on;
              end
          end
     end
end

%绘制某时刻MEO位置分布
if flag == 1
    for j = 1:1:62 
        if isempty(sate{j,t})
            continue;
        end
        if(sate_class(sate{j,t}(2),sate{j,t}(4)) == 2)
            covert_point(sate{j,t}(2)^2,sate{j,t}(1),sate{j,t}(4),sate{j,t}(3),sate{j,t}(5),sate{j,t}(17),'red');
            hold on;
        end
    end
end

        在此代码实现过程中，其输入参数flag决定了绘制MEO的类型，flag == 0，绘制所有MEO所有时刻的位置分布；flag == 1，绘制 t 时刻所有MEO的位置分布，t为输入参数。核心代码covert_point如下：

covert_point.m 

%求地球固定坐标系下轨道上一点的坐标
%输入参数：轨道六要素，颜色

function covert_point(a,e,i0,omega,w,niu,color)
omega = -omega;i0 = -i0;w = -w;
R3_omega = [cos(omega),sin(omega),0;-sin(omega),cos(omega),0;0,0,1];
R1_i0 = [1,0,0;0,cos(i0),sin(i0);0,-sin(i0),cos(i0)];
R3_w = [cos(w),sin(w),0;-sin(w),cos(w),0;0,0,1];
R = R3_w*R1_i0*R3_omega;
A = a*(1-e^2)/(1+e*cos(niu))*[cos(niu);sin(niu);0];
B = R*A;
scatter3(B(1),B(2),B(3),[],color,"pentagram");

此代码实现步骤与covert_plane相近，只不过为求地球固定坐标系下轨道的一点坐标 。

         draw_MEO_local也用到了sate_class函数，此函数作用是判断卫星的类型

sate_class.m

%判断卫星类型
%输入参数：半长轴的平方根，轨道倾角
%输出参数：sate_class == 0时为GEO
%          sate_class == 1时为IGSO
%          sate_class == 2时为MEO

function sate_class = sate_class(sqrt_a,i0)
if(sqrt_a>6000)
    if(i0<0.11)
        sate_class = 0;
    else
        sate_class = 1;
    end
else
    sate_class = 2;
end

推算两小时内igso位置和速度变化

%两小时内igso位置和速度变化
j = 1;
for i = 1:5:7200
    [local_x,local_y,local_z,velo_x,velo_y,velo_z] = draw_point_perturbation(7,1,i,"magenta");
    local(j,1) = local_x;local(j,2) = local_y;local(j,3) = local_z;
    velo(j,1) = velo_x;velo(j,2) = velo_y;velo(j,3) = velo_z;
    j = j + 1;
end

%绘制位置和速度随时间变化曲线
t = 1:5:7200;
figure
subplot(2,3,1);
plot(t,local(:,1),"red");xlabel('t');ylabel('local_x');
subplot(2,3,2);
plot(t,local(:,2),"red");xlabel('t');ylabel('local_y');
subplot(2,3,3);
plot(t,local(:,3),"red");xlabel('t');ylabel('local_z');
subplot(2,3,4);
plot(t,velo(:,1),"blue");xlabel('t');ylabel('velo_x');
subplot(2,3,5);
plot(t,velo(:,2),"blue");xlabel('t');ylabel('velo_y');
subplot(2,3,6);
plot(t,velo(:,3),"blue");xlabel('t');ylabel('velo_z');

         此部分代码通过draw_point_perturbation（）函数得到了每隔5s的位置坐标分量与速度分量。然后通过plot绘制出每个分量随时间t的变化曲线。

draw_point_perturbation.m

%绘制有摄动下地球固定坐标系下卫星位置与速度
%输入参数：flag：卫星号
%          t:卫星某一时刻
%          t_k:规划时间
%          color：颜色
%输出参数：x,y,z：卫星位置
%          a,b,c:卫星速度

function [x,y,z,a,b,c] = draw_point_perturbation(flag,t,t_k,color)
global sate;
% global local;
% global velo;
%计算信号发射时刻
% str1 = '2023/09/01 ';
% str2 = num2str(t);
% str3 = ':0:0.0';
% str = strcat(str1,str2,str3);
% delta_time = delta_StandardTime('2006/01/010:0:0.0',str);
% while 1
%     if(delta_time > 302400)
%         delta_time = delta_time - 604800;
%     else break
%     end
% end
%计算规划时间
% t_k = delta_time - sate{flag,t}(7);

%计算平均角速度
GM = 3.986004e14;
n0 = sqrt(GM/sate{flag,t}(2)^6);
n = n0 + sate{flag,t}(8);
%计算平近点角
M_k = sate{flag,t}(6) + n*t_k;
%计算偏近点角
E0 = M_k;% 初始近似解
e = sate{flag,t}(1);
f = @(E) M_k + e*sin(E);
epsilon = 1e-10; % 迭代精度要求
max_iter = 500; % 最大迭代次数
[E_k, num] = iteration_method(f, E0, epsilon, max_iter);
%计算真近点角
niu_k = 2*atan(sqrt((1+e)/(1-e))*tan(E_k/2));
%计算t时刻升交点角距
fai_k = niu_k + sate{flag,t}(5);
%计算摄动矫正项
delta1_u_k = sate{flag,t}(13)*sin(2*fai_k) + sate{flag,t}(12)*cos(2*fai_k);
delta1_r_k = sate{flag,t}(11)*sin(2*fai_k) + sate{flag,t}(16)*cos(2*fai_k);
delta1_i_k = sate{flag,t}(15)*sin(2*fai_k) + sate{flag,t}(14)*cos(2*fai_k);
%计算校正后升交点角距，轨道向径矫正项，轨道倾角矫正项
u_k = fai_k + delta1_u_k;
r_k = sate{flag,t}(2)^2 * (1 - sate{flag,t}(1) * cos(E_k)) + delta1_r_k;
i_k = sate{flag,t}(4) + sate{flag,t}(10) * t_k + delta1_i_k;
%计算升交点轨道直角坐标系坐标
x_k = r_k*cos(u_k);
y_k = r_k*sin(u_k);
%计算升交点经度与地心地固直角坐标
w_e = 7292115e-11;
%MEO  IGSO
if sate_class(sate{flag,t}(2),sate{flag,t}(4)) == 1 || sate_class(sate{flag,t}(2),sate{flag,t}(4)) == 2    
    omega_k = sate{flag,t}(3) + (sate{flag,t}(9) - w_e)*t_k - w_e*sate{flag,t}(7);
    R3__omega_k = [cos(omega_k),-sin(omega_k),0;sin(omega_k),cos(omega_k),0;0,0,1];
    R1__i_k = [1,0,0;0,cos(i_k),-sin(i_k);0,sin(i_k),cos(i_k)];
    A = R3__omega_k * R1__i_k * [x_k;y_k;0];
end
%GEO
if sate_class(sate{flag,t}(2),sate{flag,t}(4)) == 0
    omega_k = sate{flag,t}(3) + sate{flag,t}(9) * t_k - w_e*sate{flag,t}(7);
    R3__omega_k = [cos(omega_k),-sin(omega_k),0;sin(omega_k),cos(omega_k),0;0,0,1];
    R1__i_k = [1,0,0;0,cos(i_k),-sin(i_k);0,sin(i_k),cos(i_k)];
    R1__5 = [1,0,0;0,cosd(5),-sind(5);0,sind(5),cosd(5)];
    R3_w_e_t_k = [cos(w_e*t_k),sin(w_e*t_k),0;-sin(w_e*t_k),cos(w_e*t_k),0;0,0,1];
    A = R3_w_e_t_k * R1__5 * R3__omega_k * R1__i_k * [x_k;y_k;0];
end

%计算速度
M_k_dot = n;
E_k_dot = M_k_dot/(1 - e * cos(E_k));
niu_k_dot = sqrt(1-e^2)*E_k_dot/(1-e*cos(E_k));
fai_k_dot = niu_k_dot;
delta1_u_k_dot = 2*fai_k_dot*(sate{flag,t}(13)*cos(2*fai_k) - sate{flag,t}(12)*sin(2*fai_k));
delta1_r_k_dot = 2*fai_k_dot*(sate{flag,t}(11)*cos(2*fai_k) - sate{flag,t}(16)*sin(2*fai_k));
delta1_i_k_dot = 2*fai_k_dot*(sate{flag,t}(15)*cos(2*fai_k) - sate{flag,t}(14)*sin(2*fai_k));
if sate_class(sate{flag,t}(2),sate{flag,t}(4)) == 1 || sate_class(sate{flag,t}(2),sate{flag,t}(4)) == 2
    omega_k_dot = sate{flag,t}(9) - w_e;
else
    omega_k_dot = sate{flag,t}(9);
end
i_k_dot = sate{flag,t}(10) + delta1_i_k_dot;
r_k_dot = sate{flag,t}(2)^2*e*E_k_dot*sin(E_k) + delta1_r_k_dot;
u_k_dot = fai_k_dot + delta1_u_k_dot;
x_k_dot = r_k_dot*cos(u_k) - r_k*u_k_dot*sin(u_k);
y_k_dot = r_k_dot*sin(u_k) + r_k*u_k_dot*cos(u_k);
% z_k = 0;
% X_k_dot = -y_k*omega_k_dot - (y_k_dot*cos(i_k)-z_k*i_k_dot)*sin(omega_k) + x_k_dot*cos(omega_k);
% Y_k_dot = x_k_dot*sin(omega_k) + (y_k_dot*cos(i_k)-z_k*i_k_dot)*cos(omega_k) + x_k_dot*sin(omega_k);
X_k_dot = x_k_dot*cos(omega_k) - x_k*omega_k_dot*sin(omega_k)-y_k_dot*cos(i_k)*sin(omega_k) + i_k_dot*y_k*sin(i_k)*sin(omega_k)-omega_k_dot*y_k*cos(i_k)*cos(omega_k);
Y_k_dot = x_k_dot*sin(omega_k) + x_k*omega_k_dot*cos(omega_k)+y_k_dot*cos(i_k)*cos(omega_k) - i_k_dot*y_k*sin(i_k)*cos(omega_k)-omega_k_dot*y_k*cos(i_k)*sin(omega_k);
Z_k_dot = y_k_dot*sin(i_k) + y_k*i_k_dot*cos(i_k);
%MEO  IGSO
if sate_class(sate{flag,t}(2),sate{flag,t}(4)) == 1 || sate_class(sate{flag,t}(2),sate{flag,t}(4)) == 2
    B = [X_k_dot;Y_k_dot;Z_k_dot];
end
%GEO
if sate_class(sate{flag,t}(2),sate{flag,t}(4)) == 0
    R1__5 = [1,0,0;0,cosd(5),-sind(5);0,sind(5),cosd(5)];
    R3_w_e_t_k = [cos(w_e*t_k),sin(w_e*t_k),0;-sin(w_e*t_k),cos(w_e*t_k),0;0,0,1];
    B = R3_w_e_t_k * R1__5 * [X_k_dot;Y_k_dot;Z_k_dot];
end
scatter3(A(1),A(2),A(3),[],color,"diamond");
% local(t_k,1) = A(1);local(t_k,2) = A(2);local(t_k,3) = A(3);
% velo(t_k,1) = B(1);velo(t_k,2) = B(2);velo(t_k,3) = B(3);
x = A(1);y = A(2);z = A(3);
a = B(1);b = B(2);c = B(3);

        此部分代码通过输入卫星号、时刻等，自动识别卫星的种类（GEO、MEO、IGSO），然后计算卫星的位置和速度，实现方式通过转换矩阵来实现。

四、运行结果

        这里以第11号卫星——IGSO为例，我们可得运行结果如下：

        将卫星换为1号卫星——地球同步轨道卫星的第17时刻 结果如下：