当前位置：首页 > news >正文

【DeepSeek】Deepseek辅组编程-通过卫星轨道计算终端距离、相对速度和多普勒频移

news 来源：原创 2025/5/10 17:41:59

引言

笔者在前面的文章中，介绍了基于卫星轨道参数如何计算终端和卫星的距离，相对速度和多普勒频移。
【一文读懂】卫星轨道的轨道参数（六根数）和位置速度矢量转换及其在终端距离、相对速度和多普勒频移计算中的应用

Matlab程序

本文笔者将尝试通过deepseek-R1模型，编写Matlab程序，完成上文中的运算过程。
以下是根据文档内容创建的MATLAB程序，用于计算卫星与地面终端之间的距离、相对速度及多普勒频移：

%% 主计算函数
function [distance, v_rel, doppler] = satellite_calculator(...
    a, e, i, Omega, omega, nu, lat_gr, lon_gr, h_gr, t)
    % 参数校验
    validateattributes(a, {'numeric'}, {'positive'});
    % ...其他参数校验
    
    % 轨道计算
    p = a*(1 - e^2);
    r = p/(1 + e*cos(nu));
    
    % 近焦点坐标系状态
    r_pf = [r*cos(nu); r*sin(nu); 0];
    v_pf = [-sqrt(mu/p)*sin(nu); sqrt(mu/p)*(e + cos(nu)); 0];
    
    % 坐标转换
    Q = dcm_body_to_eci(omega, i, Omega);
    r_eci = Q * r_pf;
    v_eci = Q * v_pf;
    
    % ECI转ECEF
    [r_ecef_sat, v_ecef_sat] = eci_to_ecef(r_eci, v_eci, earth_omega, t);
    
    % 地面终端坐标
    [r_ecef_gr, v_ecef_gr] = geodetic_to_ecef(lat_gr, lon_gr, h_gr);
    
    % 相对参数计算
    delta_r = r_ecef_sat - r_ecef_gr;
    distance = norm(delta_r);
    
    delta_v = v_ecef_sat - v_ecef_gr;
    v_rel = dot(delta_v, delta_r)/distance;
    
    % 多普勒计算
    c = 299792458;
    f0 = 1575.42e6; % GPS L1频率
    doppler = -v_rel/c * f0;
end

function Q = dcm_body_to_eci(omega, i, Omega)  
% 生成从近焦点坐标系到ECI的旋转矩阵  
% 旋转顺序：Rz(Omega) * Rx(i) * Rz(omega)  

Rz_Omega = [cos(Omega) -sin(Omega) 0;  
            sin(Omega)  cos(Omega) 0;  
            0           0          1];  
        
Rx_i = [1       0        0;  
        0   cos(i) -sin(i);  
        0   sin(i)  cos(i)];  
        
Rz_omega = [cos(omega) -sin(omega) 0;  
           sin(omega)  cos(omega) 0;  
           0           0          1];  
       
Q = Rz_Omega * Rx_i * Rz_omega;  
end 

function [r_ecef, v_ecef] = eci_to_ecef(r_eci, v_eci, earth_omega, t)  
% 输入：  
%   t - 自参考历元的时间（秒），用于计算地球旋转角度  
% 输出：  
%   r_ecef - ECEF位置矢量 [m]  
%   v_ecef - ECEF速度矢量 [m/s]  

theta = earth_omega * t; % 地球旋转角度  

% 位置转换  
Rz_theta = [cos(theta) -sin(theta) 0;  
           sin(theta)  cos(theta) 0;  
           0           0          1];  
r_ecef = Rz_theta * r_eci;  

% 速度转换  
omega_cross = [0 -earth_omega 0;  
              earth_omega 0 0;  
              0 0 0];  
v_ecef = Rz_theta * (v_eci - omega_cross * r_eci);  
end  

function [r_ecef, v_ecef] = geodetic_to_ecef(lat, lon, h)  
% WGS-84参数  
a = 6378137.0;        % 半长轴 [m]  
f = 1/298.257223563;  % 扁率  
e2 = 2*f - f^2;       % 第一偏心率的平方  

% 地理坐标转换  
N = a / sqrt(1 - e2*sin(lat)^2);  
x = (N + h) * cos(lat) * cos(lon);  
y = (N + h) * cos(lat) * sin(lon);  
z = (N*(1 - e2) + h) * sin(lat);  

r_ecef = [x; y; z];  

% 地面终端速度（地球自转导致）  
omega_earth = 7.292115e-5;  % [rad/s]  
v_ecef = cross([0;0;omega_earth], r_ecef);  
end

功能说明

轨道参数初始化：定义卫星的轨道六要素和地球参数。轨道半径计算：根据开普勒运动公式计算当前位置的轨道半径。
位置和速度矢量：在近焦点坐标系中计算卫星的位置和速度。坐标转换：通过旋转矩阵将近焦点坐标系的位置和速度转换到ECI坐标系。
ECEF坐标转换：考虑地球自转效应，将ECI坐标转换为ECEF坐标。地面终端位置计算：将地理坐标转换为ECEF坐标。
距离、速度和多普勒频移计算：计算卫星与地面终端之间的关键参数。

测试验证

使用示例参数运行程序，验证输出结果是否合理。
检查单位是否一致（例如：距离单位为米，速度单位为米/秒）。
确保所有角度转换为弧度进行计算。

注意事项

确保所有角度参数在计算前已转换为弧度。
检查坐标系转换矩阵的正确性。
确保地球自转角速度的方向和大小与参数一致。

验证程序

以下是为验证主程序功能而创建的测试程序，包含多个测试用例和结果验证：

% 测试用例集合
test_cases = struct(
    'a', {7.0e7, 8.0e7, 7.2e7}, ...
    'e', {0.1, 0.2, 0}, ...
    'i', {45*pi/180, 60*pi/180, 30*pi/180}, ...
    'Ω', {120*pi/180, 180*pi/180, 60*pi/180}, ...
    'ω', {30*pi/180, 45*pi/180, 15*pi/180}, ...
    'ν', {0*pi/180, 90*pi/180, 180*pi/180});

% 测试结果存储
results = struct('distance', [], 'speed', [], 'doppler', []);

% 运行所有测试用例
for i = 1:length(test_cases.a)
    a = test_cases.a{i};
    e = test_cases.e{i};
    i_angle = test_cases.i{i};
    Ω = test_cases.Ω{i};
    ω = test_cases.ω{i};
    ν = test_cases.ν{i};
    
    % 调用主程序进行计算
    [distance, speed, doppler] = satellite_communication_calculator(...
        a, e, i_angle, Ω, ω, ν);
    
    % 存储结果
    results.distance(end+1) = distance;
    results.speed(end+1) = speed;
    results.doppler(end+1) = doppler;
    
    % 显示当前测试进度
    fprintf('测试用例 %d/%d 已完成\n', i, length(test_cases.a));
end

% 显示所有测试结果
disp('所有测试结果:');
disp(results);

% 可视化结果
figure;
subplot(3,1,1);
plot(results.distance);
title('距离（米）');
xlabel('测试用例');
ylabel('距离');

subplot(3,1,2);
plot(results.speed);
title('相对速度（米/秒）');
xlabel('测试用例');
ylabel('速度');

subplot(3,1,3);
plot(results.doppler);
title('多普勒频移（Hz）');
xlabel('测试用例');
ylabel('频移');