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

news2025/2/10 22:28:03

引言

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

Matlab程序

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

% 定义轨道六要素和地球参数
a = 7.0e7;         % 半长轴（单位：米）
e = 0.1;          % 偏心率
i = 45 * pi/180; % 轨道倾角（弧度）
Ω = 120 * pi/180; % 升交点赤经（弧度）
ω = 30 * pi/180; % 近地点幅角（弧度）
ν = 0 * pi/180;   % 真近点角（弧度）

μ = 3.986e14;     % 地球标准重力参数（m³/s²）
R = 6378137;      % 地球赤道半径（米）
earth_ω = 7.2921e-5; % 地球自转角速度（rad/s）

% 计算轨道半径
r = a * (1 - e * cos(ν));

% 近焦点坐标系下的位置和速度矢量
r_perifocal = [r * cos(ν), r * sin(ν), 0];
v_perifocal = [-sqrt(μ * a) * e * sin(ν) / (r), 
              sqrt(μ * a) * (1 + e * cos(ν)) / r, 
              0];

% 构建旋转矩阵
Q = dcm_body_to_eci(ω, i, Ω);

% 轨道坐标系到ECI坐标系的转换
r_eci = Q * r_perifocal';
v_eci = Q * v_perifocal';

% 转换到ECEF坐标系
[lat, lon, h] = eci_to_ecef(r_eci, v_eci, earth_ω);

% 地面终端位置（示例：北京）
lat0 = 39.9042 * pi/180; % 纬度
lon0 = 116.4074 * pi/180; % 经度
h0 = 50;                   % 高度（米）

% 地理坐标转ECEF坐标
[r_ecef, v_ecef] = geodetic_to_ecef(lat0, lon0, h0);

% 计算距离
distance = norm(r_eci - r_ecef);

% 计算相对速度
velocity_diff = v_eci - v_ecef;
relative_speed = norm(velocity_diff);

% 计算多普勒频移
c = 299792458; % 光速
f0 = 1.5e9;    % 载波频率（示例值）
doppler_shift = (relative_speed / c) * f0;

% 输出结果
fprintf('卫星与地面终端的距离：%.2f km\n', distance / 1000);
fprintf('相对速度：%.2f m/s\n', relative_speed);
fprintf('多普勒频移：%.2f Hz\n', doppler_shift);

功能说明

轨道参数初始化：定义卫星的轨道六要素和地球参数。轨道半径计算：根据开普勒运动公式计算当前位置的轨道半径。
位置和速度矢量：在近焦点坐标系中计算卫星的位置和速度。坐标转换：通过旋转矩阵将近焦点坐标系的位置和速度转换到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('频移');